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AERO POWER 


Bristol Siddeley 
lift/thrust engine chosen 


for Fokker 


supersonic VTOL aircraft 


A Bristol Siddeley lift/thrust turbofan engine 
has been chosen to power the very advanced 
Fokker/Republic VTOL aircraft with variable 
wing geometry. This aircraft will have long 
range and high supersonic all-weather capabili- 
ties to carry out its strike/reconnaissance duties. 
The engine for Fokker/ Republic is an ad- 
vanced version of the Pegasus lift/thrust turbo- 
fan which already powers the Hawker P 1127 
and is supported by the Mutual Weapons 
Development Programme for NATO. 


PEGASUS SUBSTANTIATES 
NEW LIFT/THRUST CONCEPT 
The rapid development and outstanding per- 
formance of the Pegasus turbofan has proved 
the soundness of the Bristol Siddeley lift/thrust 
philosophy. This lays down that the total in- 
stalled thrust of an aircraft's main powerplant 
must be available for vertical take-off and land- 
ing, hovering, as well as for forward flight. 

The revolutionary propulsive system of the 


Pegasus realises this basic concept and at the 
same time provides great reserves of power in 
military applications. The engine is equally 
suitable for single or multi-engined, subsonic or 
supersonic aircraft. In applications which re- 
quire very high payload or extreme range. 
Pegasus turbofans can be used in conjunction 
with pure lifting engines. 


UNIQUE DESIGN—UNIQUE PERFORMANCE 
The unique feature of the Pegasus is the use of 
four movable nozzles, evenly disposed around 
the centre of gravity, which give directional 
control to its total thrust. The engine’s high by 
pass ratio results in an excellent thrust/weight 
ratio, a low fuel consumption and a low noise 
level. Installation is conventional, and operation 
is simple and independent of all fixed ground 
equipment because the movable nozzles permit 
the problems of ground erosion, recirculation 
and ingestion of debris to be reduced to a mini 
mum. 


BRISTOL SIDDELEY ENGINES LIMITED Turbojets Turbotans Turboprops Piston Engines Ramiets 


Rocket Engines Marine and Industrial Gas Turbines 


The Aeronautical Quarterly 


Marine, Rail and Industrial Diese! Engines 


Precision Engineering Products 


November 1961 


Republic 


aeronautical problems. 


The Effects of Heat Transfer on 
Laminar-Boundary-Layer Separation 
in Supersonic Flow 


N. CURLE*, M.Sc., Ph.D. 


(Aerodynamics Division, National Physical Laboratory) 


SUMMARY : By generalising earlier work, an approximate method is derived for 
calculating compressible laminar boundary layers with an unfavourable pressure 
gradient, with the wall temperature and Prandtl number arbitrary and some 
allowance made for the viscosity-temperature law. The method merely 
requires two quadratures and uses tables of universal functions. On the basis 
of this method it is suggested that when the wall is at uniform temperature the 
viscosity-temperature law should make little difference to the results, and that, 
when the wall temperature varies spatially by a factor of two, the use of a 
linear relationship may introduce ersors of order 4 per cent in the momentum 
thickness. It is further suggested that the effects of Prandtl number differing 
from unity will increase as either the Mach number ts increased or the wall is 
cooled. These general conclusions are borne out by detailed numerical 
calculations in a few selected cases. 


The method is used to investigate “equilibrium” laminar-boundary-layer 
separation which occurs well upstream of the shock or step which induces it. 
The equations of the basic method must be solved simultaneously with a 
relationship between the (unknown) pressure gradient and the thickening of the 
boundary layer. Subject to the physically realistic approximation that the 
pressure rise to separation is a short sharp one, these equations integrate very 
easily to yield the pressure distribution in the region of interaction. It is 
predicted that the pressure coefficient, C,, at separation will equal 
0-825 (M.?—1)"'/4 R~! 4, where M, is the upstream Mach number and R is the 
Reynolds number at separation, based on distance from the leading edge and 
conditions in the main stream. Thus C, is independent of wall temperature. 
Experimental values for cases with zero heat transfer, taken from many sources, 
are scattered between 0 and 20 per cent above the theoretical value. It is also 
predicted that the pressure gradient at separation should be inversely pro- 
portional to wall temperature. The magnitude of the pressure gradient agrees 
satisfactorily with experimental values when there is no heat transfer; there 
is only limited experimental evidence concerning the magnitude of the effect of 
wall temperature, but this agrees reasonably with the predictions of the theory 


Introduction 


Interactions between shock waves and laminar boundary layers often occur in 
In many cases of practical importance the boundary layer 


*Now at the Department of Aeronautics and Astronautics, University of Southampton 
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will separate well upstream of the position where the shock strikes the layer. In 
such situations the separation is of the “equilibrium” type, in which the deflection 
of the external streamlines, caused by the boundary layer thickening, must lead to a 
pressure gradient of just the right magnitude to cause the observed thickening. 


This situation can be expressed mathematically, although the resulting analysis 
is often complicated. In principle, any method of calculating the development of a 
compressible laminar boundary layer from a knowledge of the pressure distribution 
could be used simultaneously with the relationship between thickening and pressure 
rise. Gadd: * considered just this problem by extending to compressible flow the 
relatively simple but physically well-founded method of Stratford“. An interesting 
result of Gadd’s analysis is that the predicted pressure coefficient at separation 
depends upon Mach number and Reynolds number, but not upon wall temperature 
T,,. On the other hand Gadd predicts that the extent of the region of interaction 
should be proportional to T,,°'*. The accuracy of this prediction is uncertain, as it 
is not yet firmly established how the streamwise extent of the interaction is 
dependent upon wall temperature. Gadd“ reports experiments on a flat plate in 
which there appears to be little effect of wall temperature. Later results, with 
greater heat transfer, do indicate an effect of wall temperature, and lend some 
quantitative support to the predictions of the theory although a measure of 
uncertainty remains. 


An alternative theoretical approach was suggested by Crocco and Lees’, who 
focused attention upon the mixing (or entrainment) which takes place at the edge of 
the boundary layer. Calculations of shock-wave and boundary-layer interactions 
by this method have been made by Bray, Gadd and Woodger. This work involved 
a considerable amount of computation and led to the result that cooling the wall 
reduces the streamwise extent of the region of interaction and somewhat increases 
both the pressure gradient and the pressure coefficient at the separation point. 


The present paper begins with the development of a method for calculating 
laminar boundary layers with a prescribed arbitrary adverse pressure gradient. In 
two earlier papers (Curle™ *’) such a method was presented on the assumptions that 
the Prandtl number « was unity and that the viscosity « was proportional to the 
absolute temperature T. These two limitations are relaxed in this paper. As in 
the earlier papers, an approximate integral of the thermal energy equation is 
assumed for the purpose of simplifying the solution of the momentum equation. 
The variation of » with T is assumed to be linear in a direction normal to the wall, 
but arbitrary in the streamwise direction. By this means it can be arranged that 
the viscosity is given accurately in the crucial region near to the wall, and reason- 
ably, even in the outer part of the boundary layer, where it is unimportant. The 
results take a very simple form, and the usual boundary-layer characteristics may be 
predicted by means of two quadratures, using tabulated universal functions. 


Applications to cases with prescribed pressure distributions are given in 
Section 4. Some general remarks are first made. It is a consequence of the present 
method that, for flows past a wall at uniform temperature, the viscosity-temperature 
relationship is unimportant. It is implied therefore that, in a more precise calculation, 
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it should have little effect, a result broadly in agreement with the work of 
Gadd and Young“”’. On the other hand, the use of a realistic viscosity-temperature 
law is more important for flows with non-uniform wall temperature. It is suggested 
that, if the wali temperature varies spatially by a factor of two, the errors introduced 
by a linear relationship may be of order 4 per cent as regards the momentum 
thickness and even more as regards the position of separation. It is clear that the 
linear law cannot be used when the wall temperature differs from a uniform one by 
more than this. The effects of Prandtl number « being different from unity are also 
considered. It is indicated that, for flow with zero heat transfer. separation will 
occur later when o = 0-73, say, than when « = |, a result already obtained by Gadd“. 
It is further indicated that the magnitude of this effect will increase as the wall is 
cooled or as the Mach number is increased. 


Some explicit numerical results are given for three particular problems. The 
first is for flow with linearly increasing pressure, the Mach number at the leading 
edge being 2, and the wall temperature being the free-stream temperature at the 
leading edge. The other two cases are each for flow with a linearly retarded external 
velocity, the Mach number at the leading edge being 4, and the wall temperature 
being (a) the free-stream temperature at the leading edge and (hb) the free-stream 
stagnation temperature. Solutions with o=1 and « CC T were given in Ref. 7, case 
(b) having zero heat transfer when ~=1. Since uniform wall temperatures are 
considered in each case, the assumption that »« OC T is retained. Further calcula- 
tions with «7=0-73 are made, however, and are compared with the «=I results, 
showing changes in separation distance of order 3 per cent, 7 per cent and 0 per cent 
respectively due to the change in Prandtl number, with appropriate changes in the 
main boundary-layer characteristics. These results bear out the conclusion already 
mentioned that the effect becomes more marked as the Mach number increases and 
as the wall ts cooled. 


Application is then made to the shock-wave and boundary-layer interaction 
problem. It is found experimentally that, when the upstream Mach number M, and 
the Reynolds number at separation, R=wu,x/v,, are in the ranges 2=M, <4. 
10° =< R<10°, the value of u, at separation is only 1-2 per cent less than its 
upstream value. Accordingly, the assumption is made that the pressure rise to 
separation is a short sharp one, so that du,/dx is non-zero but u, is nearly constant. 
With no further approximations, the equations of Refs. 7 and 8, together with an 
equation relating the pressure rise to the thickening of the boundary layer, are 
reduced to a single equation for the pressure as a function of x, and this equation is 
integrated very easily. The solution can be continued beyond separation, provided 
that the relevant universal functions in Ref. 7 are defined in this region; this is 
considered in Appendix III. 


The method predicts, with the earlier work of Gadd, that the pressure 
coefficient at separation is independent of wall temperature. It also predicts that 
the extent of the region of interaction is directly proportional to T., so that the 
pressure gradient at separation is inversely proportional to T,. Quantitatively, the 
predicted pressure coefficient at separation is C,=0°825 (M,?—-1)''* R-''*. This 
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agrees well with experimental results for zero heat transfer, which are scattered from 
0 to 20 per cent above this value, and compares favourably with other theoretical 
predictions. The accuracy of the prediction concerning the pressure gradient at 
separation is, as implied earlier, still uncertain, although the results of the most 
recent experiments are encouraging in this respect. 

Numerical results have been obtained in a few cases for which experimental 
results were available. For a case of two-dimensional flow with zero heat transfer, 
with M,=2 and R=7~ 10‘, the predicted pressure distribution agrees well with 
experiment in the region upstream of where the shock strikes the boundary layer. 
Some published results were also available (Gadd) for three cases with M,=3 and 
R=4:2~ 10°. In one case there was zero heat transfer, and in the others the wall 
was cooled to a temperature 7, =0°88 T, and heated to T,,=1:25 T,. The predicted 
pressure distribution agrees remarkably well with experiment for the case of zero 
heat transfer. The theory does, however. predict an effect of wall temperature which, 
though small, is greater than that found in the experiments. Later experiments, 
with M,=2-70 and R=1-:5~x 10°, with zero heat transfer in one case and the wall 
heated to 7, = 1-5 T, in another, yield results which lend some support to the theory 
so that, although no firm conclusions can be reached as to the accuracy of the 
predicted effect of wall temperature, the theory is at least qualitatively correct in 
this respect. 


NOTATION 

x,y distances parallel and normal to the wall y=0 

u,v velocity components parallel to x and y axes, respectively 
p pressure 
p density 
T absolute temperature 
u Viscosity 
v kinematic viscosity, 
a_ speed of sound 
M Mach number, =u/a 
o0 Prandtl number 

6,,6, boundary layer displacement and momentum thickness, respectively 

6,’,%,",6,’ defined in equations (25), (26), (29) 

R_ Reynolds number, =u,x/», 


C, specific heat at constant pressure 


C, pressure coefficient, 
\ 

T, total temperature, — T + 
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K_ defined in equation (3) 
6, defined in equation (5) 


€=1—u,7/u 


1 
Uv stream function, pu dy 
Ms. 


Y related to y as in equation (16) 

\ felated to ¥ as in equation (17) 

Ca function of x, appearing in equation (19) 
G_ defined in equation (23) 


Q defined in equation (31) 


H’ =3,'/3,' 
ay) 

m Te 7 
L’=2m’ (H’ + 2)4+2l 

G, =H’ —2:61 


c areference length 
t defined by equation (53) 

F.X,A defined by equation (61) 

Suffixes 

0 conditions at edge of boundary layer upstream of interaction* 
1 local conditions at edge of boundary layer 
s standard reference conditions 
w conditions at the wall 


Zz conditions appropriate to zero heat transfer 


2. Temperature Profile 


To make it possible to perform a relatively straightforward integration of the 
momentum equation it is necessary first to make some assumption or approximation 
regarding the temperature profile. This is now considered. 


It is well known that provided (i) the Prandtl number « is unity, (ii) there is no 
pressure gradient, (iii) the wall temperature 7, is constant, then the temperature 
profile takes the form 


Ta T+ T.) 


l 
*Except in F. defined by equation (65) 
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This will also be the profile locally at a sharp leading edge, whatever the pressure 
gradient and wall-temperature distribution. It is extremely difficult, except in the 
special case for which equation (1) is valid, to calculate the downstream development 
of the temperature profile, for which reason the author”? suggested a method of 
choosing a quadratic in u/u, as an approximation to the total temperature Ty. 
This approximation takes the form 


u 


Tu =Te+(T,-Tw) +KT, (“ - (2) 
u 

T,-—Ty(0 
where KT, =T,, (0)-T.+ (3) 

1+6,(€-1) 
with g=1- (4) 

T,, (0)/T, -—1+4(y¥-1) Me? 


T..(0)/T,+1+4(y-1) M2? 


Throughout this analysis it is assumed that the boundary layer starts at a sharp 
leading edge, at which u,=u,. For flow past a bluff body, with the boundary layer 
developing from a forward stagnation point, the analysis may be applied down- 
stream of the pressure minimum provided that the origin x=0 is taken there 


rhe use of this profile leads to a marked improvement in calculations of the 
velocity boundary layer in an adverse pressure gradient, as compared with calcula- 
tions assuming that the temperature profile takes the local flat-plate value (equation 
(1)). In typical cases the error in predicting the position of separation appears to 
be halved. 

When the Prandtl number is not equal to unity, the position is further 
complicated and it is necessary to generalise equation (2) further, for which reason 
we require first a generalisation of the flat-plate profile (equation (1)). A major step 


in this direction was made by Monaghan‘''’, who wrote equation (1) as 
T =T.+(T,—Tw) ——-M,?T, (6) 
2 
and showed that a similar form 
ui, 2 


yields an accurate prediction of the temperature field for values of u/u, up to about 
0-8 when o =0:725 and «=1:25. To obtain a formula for T which also satisfies 
the correct boundary condition at the edge of the boundary layer, we shal! add a 
term A (x) u*/u,* to the right-hand side of equation (7) and choose A (x) so that 
T—>T,asu—>u,. This gives the result 


T=T,+(T,-T,) (1 +o"! (— - - 
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We must now add a term, analogous to the last term in equation (2) to account for 
downstream variations in equation (8). In doing this it is reasonable to take this 
term, KT, [(u/u,)—(u*/u,*)], multiply the two terms in the bracket by o'/* and « 
respectively and add on a suitable multiple of u’/u,* to satisfy the external boundary 
condition. This procedure yields 
T -T,~(T.-T,) (1— —,) (T.-T.+KT,) - 

u, 
u (9) 


3 


T, (K+ (_, 
2 ou 
Finally a minor amendment must be made in the definition of K if incon- 
sistencies are to be avoided. From equation (9) it is seen that K must be zero when 
there is no heat transfer, which occurs when the wall temperature distribution is 


| 2 

The definition of K given by equation (3) is easily seen to satisfy this condition only 
when the wall temperature for zero heat transfer is a constant, that is when o =1. 
For general values of « we add to equation (3) a function B(é,), such that K =0 
when 7, =T This leads to the following result, accepted in all that follows, 
namely 


KT, ‘ (11) 


1+4,(€-1)) 

It should perhaps be remarked that the temperature profile given by equations (9) 
and (11) is at best empirical, and can be finally justified only by the conclusions 
drawn from it. In this respect the agreement with calculations by Gadd‘ and 
Young”, discussed later, is encouraging. 


3. Transformation of the Momentum Equation 


For steady two-dimensional compressible flow the boundary-layer momentum 
equations and the equation of continuity take the form 


Ou Ou 1 o Ou 
Cx Cy x p CV) 
Op 
(13) 
fal 
(pu)+ ~- (pv)=0. (14) 
Xx oy 


Equation (14) is satisfied by introduction of a stream function v, with 
uU= Pp, » = (15) 
dy f ax 

suffix s denoting some arbitrary reference condition. 
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Following Howarth”, we change from independent variables x, v to x, Y, 


y 


[Tas 
where Y= dy, 
0 


Then, exactly as in Howarth’s analysis, it is quite straightforward to show that 


CU dx 07x dx 1 dp 07y (18) 
Ox Oy Ox0oYOY OY* Ox 2pdcx 


To deai with the viscous term in (12), it is assumed that the viscosity is related 
to the temperature by the law 


Ctx) 


Thus, for fixed x, » varies linearly with 7 across the boundary layer, but an 
arbitrary variation in the streamwise direction is permitted. Following the suggestion 
of Chapman and Rubesin®, C (x) may be fixed so that the viscosity is given either 
exactly at the wall or fairly accurately throughout the region near to the wall where 
viscosity is most important. 


Accepting equation (19), and leaving C (x) to be determined in one of these 
ways, it then follows that 


T. 


(20) 
Ps oY’ 
p | oy Oy) Dp, 
wo 
(21) 
Gy or” 


Then, substituting from equations (18) and (21) into equation (12), we obtain 


0*y Cy O7y Oy Op l Cp 
OxOYOY 0x 2pdox” p ax 
lu x 
following Howarth, with G (23) 
316 
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We now integrate across the boundary layer, from Y=0 to Y=00. Upon 
defining 


~ 


| (1- ay (24) 
u, 
and | (1- dY, (25) 
J uy, 
this yields 
+ 26, )=Cv, | oY). Mi ay (G 1) dY. (26) 
From equation (23) we deduce that 
(G—1)dY ) d¥ — 4,M,,73,’, (27) 
and, upon substituting from equation (9) for T, this yields 
(G 1) dy ( +- j+o 0, — 
, —] , , 
where | (1 dY (29) 
u, 


By substituting from equation (28) into (26), and writing A, for the ratio 6,’/4,’, we 
find that 


 du,|{Tw,,, 
ax \T, 


Cy, (30) 


where Q=2-4M,7+ h, (l-o'*)+K {hyo' (h,—1)o} + 


It is noted that Q takes the form of the sum of a constant plus multiples of 
M,’, ((T./T,)—1] and K, the multiples being respectively — 4, 0, 1, when o is unity. 
When o takes other values it is necessary to know h, before Q can be tabulated. It 
was pointed out by the author, following Truckenbrodt“”, that A, is almost 
constant, and that a reasonable mean value is 1:58. It is worth noting that, with 
this value of h,, the coefficient of K in (31) becomes 1-58 o'/* —0-58 o, which lies 
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TABLE I 
THE FUNCTIONS H’(m’) AND ['(m’) 

m H’ m H’ 
2-00 0: 500 0-060 2:99 0-113 0.075 4°65 
0°20 2:07 0-463 0°064 3-04 0-104 0-070 484 -0'116 
2°18 0404 0:068 3:09 0-095 0-065 5:05 0-118 
0:382 0-072 0-085 0-060 5°30 
2:28 0-359 0-076 3°22 0072 0-055 5-60 0-115 
0-080 2:34 0:333 0-080 3-30 0-056 0-050 5:96 0-112 
0-064 2-39 0-313 0-084 3°39 0-038 0-045 6°41 0-108 
2:44 0-291 0-086 344 0-027 0-040 7-00 —0:103 
0-032 2:49 0 268 0-O88 3-49 O-O1S 0-035 7°74 0-097 
O-OlF 2°55 (244 0-090 3°55 0 0-030 8°64 -0°089 
0 2 61 (220 0 0:025 967 -0-O8 
0-016 2°67 0-195 0-088 3-92 —()-050 0-020 10°92 —():072 
275 0°168 0-086 412 -0:070 0-015 12°61 -0:062 
0-040 281 0-153 0-084 4:26 —0:084 0-010 15:26 —0:050 
2°87 0-082 4:38 —0:094 0-005 21:05 —0°034 

0-056 2:94 0°122 0-080 4:48 0 x 0 


between 0:990 and 1-004 for values of « in the range 0:'636 = o« = 1-135, and equals 
0-999 when « =0:73, close to the value for air. For this reason it would appear 
worthwhile simply to replace the coefficient of K in (31) by unity throughout this 
range of o. 


We now define H’ . 
(33) 
NOY 
ar m Cy, dx 


so that equation (30) may be written as 


= LL’ —4m’+2m’Q ; . (35) 
where L’= 2m’ (H’ +2)+2YI. (36) 


It remains to prescribe values for the functions I’, H’, L’, in terms of m’ and 
such other variables as may be relevant, and then to transform from the artificially 
introduced 4,’, 4,’, Y, to the physically more important 5,, 6,, y. As the analysis is 
rather similar to that for the case e=1, » OC T, given in Ref. 7, these results are 
derived in Appendices I and II and the results are now used. 


It is indicated in Appendix I that H’, I’ and L’ depend primarily upon m’, being 
functionally independent of wall temperature, Mach number and Prandtl number. 
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It follows that, in general, the best results should be given by using Thwaites’s 
values, which are quoted for convenience in Table I, as modified by Curle and 
Skan©”. 


Now Thwaites showed that L’ is given, to a very good approximation, as 


L’=0°45 + 6m’, (37) 
so that equation (35) becomes 
u, [6,7 T, 
0-4: 2m’ 7). (38) 
which may be written as 
d 2u,’ (T 6,7 0:45C 
4 = + (39 


The integrating factor for this equation is 


x 


+O d. ‘ (40) 
G,=exp | | T. x 
0:45¥. G.C 
G, ul 


In Appendix II it is shown that 


6,73 (42) 
so equation (41) may be written 


0:45y, G.C 
5,7 : ; (43) 
“Gola 


It follows similarly from equations (34) and (42) that 


m=—-— u,’ (44) 


It is also shown in Appendix II that 


é, H'+ -1) yh,+K+—5— M,? [hy -(h, (45) 


the coefficient of K being almost exactly unity when 0-636 < oc < 1-135. 


Ou Ou\ (oY 
Finally, (59). (Ss). 
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After some algebra equation (46) may be reduced to 


where the last bracket will be unity if C is chosen to give the viscosity exactly at 
the wall. 


The solution of the compressible laminar-boundary-layer equations is then 
given by equations (43), (44), (45) and (47), with equations (11), (31) and (40) being 
used to determine the function G,. It will be noted that, even with the added 
generalisations of this paper, the method requires only two simple quadratures and 
the use of tabulated universal functions. 


4. Applications to Some Flows with Prescribed Pressure Gradient 


We begin with some general observations as to the effects of Prandtl number 
and viscosity-temperature relationship on the predicted laminar-boundary-layer 
development. 


Consider first the effects of the viscosity-temperature relationship, which enters 
through the function C (x), 7 being assumed to equal unity. It seems reasonable to 
suppose, following Chapman and Rubesin’s remarks, that it is most important to 
ensure that viscosity is accurately represented in the region close to the wall. If we 
go all the way and determine C (x) by the condition that viscosity is given exactly 
at the wall itself, then we find 


po (x) T, 


) 


It follows immediately that, when the wall temperature is uniform, C is a constant, 
and the results obtained by the present method will be identical with those for 
uOCT. We conclude that in a more precise calculation the viscosity-temperature 
relationship will have little effect for flows past a wall at uniform temperature. 


his conclusion agrees well with some results of Gadd“, who considers flow 
past an insulated wall with a uniform adverse external velocity gradient, and assumes 
u OCT’. He finds that the distance to separation when » =8/9 is 2 per cent greater 
than when w=1, the Mach number being M,=4. Similarly, when M,=10 the 
separation distance is 6 per cent greater when » =8/9. The corresponding results 
given by the method of Young” are also quoted by Gadd, and are 2 per cent and 
3 per cent respectively. The effect of changing » is small, even at high Mach 
number, which lends some support to the use of the approximate viscosity- 
temperature relationship of this paper. 


The position is somewhat more complicated when the wall temperature is a 
function of position, but a rough idea as to the errors introduced in taking « Oc T 
may be obtained as follows. With » oc T° equation (48) may be written 

Cwatt. (49) 
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Hence, comparing the cases »=1, » =8/9, and remembering that the integrand in 
equation (43) is always positive, we note that there is an additional factor T,~'/° 
in the integrand when » =8/9, associated with a factor T,'/* outside the integral. 
If 7. (x) increases with x, the integral will be multiplied by a factor lying between 
{T~(O)}~*/° and {7,, (x)}~'*, so the value of 6,? (and hence of m’) is multiplied by 
a factor lying between unity and {7, (x)/T, (0)}'°. Since separation occurs at a 
prescribed value of m’, the distance to separation will be decreased, by a factor of 
order (x)/T, at most. Accordingly, when the wall temperature doubles 
or halves itself, the use of a linear viscosity-temperature law may affect the predicted 
momentum thickness by about 4 per cent, and the separation position by up to 
8 per cent, with correspondingly smaller differences in the upstream characteristics 
of the boundary layer. 

Consider now the effects of Prandtl number, assuming «Oc T. These effects 


enter through the function Q (x) used in the preliminary quadrature, and defined in 
equation (31). For values of « near to unity it is easy to expand Q (x) in the forin 


(50) 
where QO, (x) - 0°53 (= 1) -0-74 
5 l 
0:53 0-21 M,’. (51) 


Unless T, is greater than T, by a factor which increases from unity at low Mach 
numbers to 1-4 at high Mach numbers, it follows that Q, will be negative. Unless 
T., is sufficiently great, G, will thus be larger when « =0-73 than when o=1, since 
retarded flow with u,’ <0 is being considered. It is deduced from equations (43) 
and (44) that m’ (x) is usually smaller when « =0-°73, since the increase in G, will 
change the integral proportionately less than the maximum change in G,. It follows 
that separation is delayed. We deduce that the effect should be considerably greater 
for a cooled wall and at higher Mach numbers. 


These general tendencies agree with the predictions of Gadd. At M,=4, for 
the linearly retarded external velocity and no heat transfer, the distance to separa- 
tion in the case 7 =0-73 is increased (relative to the case 7 = 1) by about 7 per cent, 
according to both Gadd’s method and that of Young. At M,=10 the effect is 
increased to 11 per cent according to Young’s method or 24 per cent according to 
Gadd’s method. 


Numerical results for three cases, obtained by the present method, are now 
considered. The first case is of flow with a linearly increasing pressure 


p=p, (1+*). . 


the Mach number at the leading edge being 2, and the wall being cooled to the 
temperature 7, of the main stream at the leading edge. Results of calculations with 
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FiGuRE 2. Skin friction for the case T, =T,, u. =u, (1—x/c) at Mach number 4 
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TABLE Il TABLE Ill 
SOLUTION FOR SOLUTION FOR 
ty 
vic 
0-328 0 0-328 
0:316 O-31¢€ 0-025 0404 ():407 
0:302 0°303 OOS 0475 0:483 
(3 0-288 0-289 0-075 §38 0°55] 
04 0:272 0-274 Ol ) 0 608 
Ors 0-255 0-258 0-125 0-622 0-651 
0-6 0-237 0-241 0-15 0) 639 0-677 
0-7 O-218 0-222 Q:175 
Ors 0-198 0-203 0-2 0613 0-671 
0-175 0-182 0-225 0-569 0-638 
0-150 0-158 0-25 0-493 0-582 
0-122 0-132 0-275 0-389 0-495 
0-079 0-095 03 0-180 0-368 
13 0-018 0°044 O-311 0 
1:36 0 0-331 0 


o=1 and «=0-73 are given in Table II and shown in Fig. 1, in terms of a skin- 
friction coefficient t,, defined as 


uo 

It is noted that ft, (x) is only slightly affected by the change in «, so that the position 
of separation changes from x=1:32c when «=1 to x=1-36c when o=0-73, an 
increase of 3 per cent. If separation positions are calculated by the empirically 
improved criterion indicated in Ref. 7, so that the value of m’ at separation is a 
function of M,, then the respective positions are x =1:24c and x=1-:28c, again 3 per 
cent greater. 


In the second case the external velocity decreases linearly with distance, 


u,=u (1- =). . (54) 
the wall temperature again being 7, =7,, but at a Mach number of 4. The results 
are given in Table III and Fig. 2. The change of Prandtl number produces, in this 
case, a change in separation position from x=0-3lle to x=0-33le or (with the 
empirical correction) from x=0:248c to 0:266c, an increase of about 7 per cent. 
The detailed distribution of skin friction is changed appropriately. 

In the third case the external velocity is as in the second case, namely equation 
(54), with a leading edge Mach number of 4, but the wall is at the free-stream 
stagnation temperature. It follows that there is zero heat transfer when o=1, the 
wall being slightly heated when ~=0-73. It can be seen from Table IV that the 
distribution of skin friction and the position of separation are virtually unaffected 
323 
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TABLE IV 
y-1 : 
SOLUTION FOR T,=T (1 M,?}, u, =u, (1- ,M,=4 
2 
x —— 
7=0°73 

0:328* 0°328 
0-01 0-318 0318 
0-02 0°301 0°302 *These two values were 
0°03 0:277 0-278 incorrectly given in Ref. 7. 
0:04 0:244 
0-05 0:199 
0-06 0:112* 
0:0670 - 0 
0:0671 0 


by the change in«. The reason for this is that a reduction in 7 would tend to cause 
a reduced momentum thickness and a later separation if there were zero heat 
transfer throughout. On the other hand, though there is zero heat transfer when 
o=1, the wall is slightly heated above the zero heat transfer value when o> < 1. 
This causes an opposite effect, and the two appear to balance. 


A comparison of the results obtained in these three cases bears out the general 
conclusions already mentioned, that according to the present method the effects of 
change of Prandtl number become more marked as the Mach number is increased, 
or as the wall is cooled below the zero-heat-transfer value. 


5. Application to Interactions between Shock Waves and Boundary 
Layers 

The type of situation to be considered here is that illustrated in Figs. 3 and 4, 
in which laminar-boundary-layer separation in supersonic flow is provoked by an 
oblique shock wave, a wedge on the wall, or some similar mechanism. In such a 
flow separation will usually occur well upstream of the agency which produces it, 
and will be of the “equilibrium” type, in which the deflection of the external 
streamlines due to the thickening of the boundary layer produces a pressure gradient 
which must be just sufficient to cause the requisite amount of thickening. 


CT 
DIREC 


FiGURE 3. Boundary layer 
and externally generated 
shock 


BOUNDARY LAYER 
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COMPRESSION 
WAVES 


Ficure 4. Shock generated 
from within boundary layer 
by wedge on surface 


_FLOW 
~—BIRECTION 


BOUNDARY LAYER 


Mathematically there are two relationships between the thickening of the 
boundary layer and the pressure distribution, which must be such that these relation- 
ships are consistent. On the one hand, if the pressure p(x) or the free-stream 
velocity u,(x) were known, then the usual boundary-layer characteristics could be 
calculated by the theory of Section 3, in particular 6, (x), the displacement thickness. 
On the other hand the pressure rise is related to the angle through which the 
external streamlines are deflected, that is to 

(x) (x)! +226, . . 
dx | } 
where 4,, (x) is the displacement thickness when there is no pressure gradient and 
the temperature boundary conditions are unaltered, and <=0 upstream of the shock 
or wedge apex and <=1 downstream thereof. This relationship, expressed 


mathematically, is 


5, (x) (x) + 220 (x-x,) 
/ \P , 

where x=.x, is the shock or wedge apex. Simultaneous solution of (56) and the 

various equations of Section 3 is, in general, a matter for high-speed computing 

machines. With suitable restrictions, however, an analytic solution can be 

obtained, as follows. 


Consider flows in which the pressure rise to separation is a sharp one, occurring 
over only a short distance. In such circumstances it is permissible to regard 
quantities such as the pressure, external velocity and temperature, and wall 
temperature, as being sensibly constant, although their gradients may not be set 
equal to zero. From a practical point of view this approximation seems good, for 
if the Mach number M, and the Reynolds number at separation R are in the ranges 
2<M,=<4, 10°< R =< 10*, the value of u, at separation differs by only | or 2 per 
cent from its upstream value, and in a typical case might decrease by only a further 
1 or 2 per cent between separation and the point where the shock strikes the 
boundary layer. Accordingly the value at separation is a good mean value to take. 
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With this approximation it is found that the relevant equations for the 
development of the boundary layer simplify as follows. In the region upstream of 
the interaction, u,’ is zero and hence, by (40), G, is unity. In the region of inter- 
action, (T,/7,)+@Q remains sensibly constant (provided that T, does not vary too 
rapidly with x), so that G, becomes equal to a power of u,/u,, and therefore remains 
close to unity. If, following the arguments advanced in Section 4, C is taken to be 
constant, then it can be deduced from (43) that 6,* is almost constant through the 
interaction and may be taken as 


(57) 
It then follows from (44) that 
, x , ) 
m’ (x)= —0°45 u’() E58) 
and then, from (45), 
(59) 


where all values except H’ are to be taken at separation. In the same way we may 
deduce that the displacement thickness 4,,, in the absence of a pressure gradient, is 
given by the same expression with H’ (m’) replaced by H’ (0). Accordingly it follows 
that (56) becomes 


x 


1) ds (0-45 


(M,? —1)'/? v 


*) | { H’ (m’) — 2-61} (x —x,). 


1 


(60) 


It is worth noting that this relationship is independent of «, the Prandtl number 
The reason for this is that, although (59) shows that the displacement thickness 4 
is dependent on «, the increment in 6, due to the pressure gradient is independent 
of it. 


We may simplify (60), and isolate the relevant parameters, by writing 


X=1-820 (M,?—1)'4 (=) (+ -1) (6) 


\ 1-638 6 (M,?-1)-"/4 


where R is the Reynolds number, u,x/v,, at separation. Then equations (60) and 
(61) become 
x 


F dX =H’ . (62) 
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In equation (62) m=dF/dX, . ; . (63) 


and X, is the value of X where x=x,. The lower limit of integration represents the 
leading edge x= 0, and, according to the assumptions made, is many times the inter- 
action length away, and may be taken as minus infinity. So if we define 
G, =H’ —2°61, and, if F’ denotes derivatives with respect to X, then (63) becomes 


x 
F dX =G, (F’)+A(X-X,). . (64) 


A solution of this equation may be obtained by writing 


F(X)=F,(X)+ Ag, (65) 
x 
whence F,, satisfies F, ax G, (F,’). . (66) 
Upon multiplying by F,”, and integrating from X 00 to X, we find 
Fy x x x 
x 
F,’ | F, dX —-4F,?. . (67) 
x 
Substitution for F,, dX from equation (66) gives 
F, 
F,?=2F,'G, (F,’)—2 | G,(z)dz=P(F,’) (68) 
o . (69) 


For each of these two solutions it is a simple matter to derive F,’ in terms of F,, by 


interpolation, and then to relate F, and X by 


(70) 


In proceeding beyond separation it is necessary to define the function H’ (m’) for 
such a region. Details of the best choice are given in Appendix III and Table I 
Accepting these values, the integration of (70) leads to the results given in Table V. 
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TABLE V 
THE FUNCTION F(X) OCCURRING IN SHOCK-WAVE AND BOUNDARY-LAYER INTERACTIONS 


X F(X) F’ (X) Xx F (X) F’ (xX) 
— 7-03 0:02 00103 0°70 0-40 00885 
ak 0:05 0:0237 1:86 0°50 0:0828 
—409 0-08 00351 3°21 0-60 0°0645 

3°14 0-12 5:03 0:70 00465 
0:17 00612 7°61 0-80 00323 
0:23 00736 0-90 00174 
0:29 0:0832 15°32 0-95 0-0101 

0 0:0900 23°33 1-00 0:0042 

x 1:03 0 


Since the pressure rises to separation it follows that F must be positive. 
Upstream of separation, where ¢=0, it follows that F=F,, so the positive sign in 
(69) must be taken. Thus, at separation, where F’ (Y)=0-090, the value of F is 


Fug 
It follows from (61), remembering that C, ~ 1 —(u,?/u,?), that 
C, (M,?-1)'/* (72) 


so that the pressure rise to separation depends upon Mach number and Reynolds 
number, but not upon wall temperature. 


At the position where the shock strikes the boundary layer or at the wedge 
apex, the pressure is continuous, so F is continuous. Hence, since A« is discon- 
tinuous at this point, it follows that the opposite sign in (69) must be taken for F, 
downstream of separation. This can only be so if F, changes from +4A to —43A. 
It follows that the pressure is symmetrical about x=.x,, and that the pressure rise 
to this point is one half of the final value. 


The predicted pressure coefficient at separation, equation (72), may be 
compared with experiment. Chapman, Kuehn and Larson” have correlated a 
considerable number of experimental results at zero heat transfer over a range of 
Mach numbers 1:1 <M, <3°6 and Reynolds numbers 10* << R < 10°, and find 
that C,(M,’—1)''* R''* has a mean value at separation of about 0-93, the scatter 
being of order 10 per cent. Some of the theoretical values known to the author are 
1-05 and 1:15 (Hakkinen ef al.°”), 1:03 (Gadd), 1:13 (Gadd®’) and 0-94 (Bray 
et al.). 


On the other hand the scale in the x-direction, which determines the upstream 
influence, is predicted by (61) to be directly proportional to T,, so that the pressure 
gradient at separation is inversely proportional to T,. Now some experimental 
results of Gadd“ suggest that the pressure distribution is dependent upon wall 
temperature to a remarkably small degree. On the other hand, later experiments 
(as yet unpublished) have shown that there is some dependence provided that heat 
transfer rates are sufficiently high. 
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Direct comparisons with detailed experimental pressure distributions have been 
made for a few cases. Firstly we consider some unpublished results of Gadd for 
flow along a flat plate with an upstream Mach number M, =2, and with zero heat 
transfer at the wall. By matching the theoretical and experimental pressures at the 
theoretical separation point, it is found that separation occurs 0°78 in. from the 
leading edge, with R=6-93 = 10*. It then follows from (61) that 


p/p, = 1+0°320F | 


(73) 
x= 0:78 (1+0-0463X) | 


with x measured in inches, whence Table V yields the pressure distribution. The 
results are given in Fig. 5. The agreement is good as far as the position where the 
shock strikes the boundary layer. Downstream of this position the agreement is less 
good, the experimental pressure distribution having a more pronounced foot. 


We now consider some results of Gadd, indicating the effects of heat transfer. 
These experiments were carried out on a flat plate with M,=3 and R=4-2~ 10°, 
with zero heat transfer in one case and the wall heated to a temperature 
T~=1-25 T, and cooled to T,, = 0°88 T, in the other cases. Using the same matching 
technique as before, the theoretical pressure distributions were predicted and are 
shown in Fig. 6 with the experimental values. The agreement is excellent for zero 
heat transfer. The predicted effects of heat transfer are small, but are rather greater 
than given by experiment. 


The amount by which the wall temperature differed from the zero-heat-transfer 
value was quite small in the cases just considered. Accordingly, comparison has 
also been made between the theory and some later results (as of Gadd and 
Attridge"”). These relate to flow with an upstream Mach number M’= 2-70 and 
with R=1-5 x 10°, the wall temperature being 7, =T7, in one case and T,=1°'5 T, 
in the other. The results are given in Fig. 7. theory and experiment being matched 


2b 
THEORY 
Pr © O00 EXPERIMENT 
xiO°PER INCH 
A=POSITION OF SEPARATION 
B=POSITION WHERE SHOCK STRIKES LAYER 
C=POSITION OF RE-ATTACHMENT 
ob 


Ficure 5. Theoretical and experimental pressure distribution for M,=2, R=7X 104. 
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FiGurE 6. Experimental and theoretical results for M, =3, R=42 X 10°. 
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ZERO HEAT TRANSFER T= T, 
| © HEATED WALL 
P/Py | THEORETICAL RESULTS 
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FiGuURE 7. Experimental and theoretical results for M.=2-7, R=1:5X 10°. 


at the apex of the wedge which provoked separation. We note that the agreement 
is reasonable in that the upstream influence is greater when the wall is hotter, and 
this lends valuable support to the theory; although a detailed comparison of the 
experimental and theoretical pressure distributions is a little disappointing, it should 
be remembered that, on the one hand the accuracy of the experimental pressures is 
probably not high and, on the other hand, the overall pressure rise of 60 per cent 
is rather large. Further experiments with a lower pressure rise and more accurate 
measurements would be desirable. 
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Appendix I 
FUNCTIONAL RELATIONSHIPS 
It was pointed out by the author that, when +~=1 and uc T, the functions H’, I’ 
and L’ depend mainly on m’, and are essentially independent of wall temperature and 
Mach number. It is the purpose here to examine whether a similar result holds for 
arbitrary fixed « and general viscosity-temperature law. 


It is assumed that the velocity profile may be approximated by a suitable function 
of the (transformed) normal co-ordinate Y. Thus 


=f (»), n= (74) 


— 


where the edge of the boundary layer is taken at »=1 if a method analogous to that 
of Pohlhausen''®’ is used, and at n=0o if a method analogous to that of Timman'' 
is used. Upon choosing f to satisfy certain of the boundary conditions relevant to 
(22), we find that at the edge of the boundary layer 


f—> |; Gt ree > 0 (75) 
At the wall, the primary boundary condition, u=0, yields 
p 


(0)=0. P ; (76) 


while, upon putting Y =0 in (22), we obtain 


du 72. 
(QO) (5) m (77) 


5 


The next boundary condition at the wall is obtained by differentiating (22) with respect 
to Y and then putting Y=0. This yields 


3 (T,—T, +KT, 
“Vax \ay dx lay). 


by (23) and (9), and this may be written 


—Te+KT, 


(0)=A’ 
f” ( \ T. 


As in Ref. 7, it is argued that the right-hand side of (78) is zero at a sharp leading 
edge (A’=0), at separation (f=0) and when there is zero heat transfer. Accordingly 
(78) is approximated by the simpler relationship 


November 196! 331 


7 
. 


N. CURLE 


The equations (75), (76), (77), (79) are formally the same as for incompressible flow, 
so that H’(m’), I'(m’) and L’(m’) are functionally independent of Mach number, 
wall temperature, Prandtl number and viscosity law, the effects of these quantities 
appearing only implicitly in the transformation of the normal co-ordinate and in the 
definitions of the functions, H’, I’ and L’. 


Appendix II 


THE TRUE BOUNDARY-LAYER THICKNESSES 


By definition $,= | (i- dy. (80) 
4, \ u, 
0 
By equation (16), dy= — 
f p 
v, 1/2 
) ay, (81) 
1/2 
so that (80) yields ) (82) 
Similarly 6,= (1- = dy 
PY, 
| + at 
oo 
= (=) | 1) (83) 


By substituting from equation (9), and carrying 
be easily shown to lead to 


out the relevant integrations, this may 


+K (o +o) 


H’+ ) h, + (5 /*h, —o fh, 1) + 


2 
using the fact that the coefficient of K is almost exactly unity when 0-636 < o < 1-135. 
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Appendix Ill 


FLOW DOWNSTREAM OF SEPARATION 


The universal functions used in Thwaites’s method were obtained from a considered 
analysis of all the exact and approximate solutions of the laminar-boundary-layer 
equations then known. To extend these functions into the region downstream of the 
separation point requires a similar analysis of the relevant solutions. 


As far as the approximate solutions are concerned, it is a fairly straightforward 
matter to calculate the consequences of formally applying them downstream of 
separation, and the functions /(m), H(m), L(m), are given in Figs. 8 to 10 for the 
methods of Pohlhausen"*® and Timman°’, and for the exact similarity solutions given 
originally for the upstream flow by Falkner and Skan‘“!) and extended to flow beyond 
separation by Stewartson'?). 


In spite of the wide divergences between the solutions, there are points of similarity. 


(i) In each case the value of m increases to a maximum, and then decreases to 
zero far iownstream of separation, corresponding physically to a considerable “ dead- 
air region. 


(ii) For the exact “ similarity * solutions the maximum of m occurs at separation, 
where there is a cusp in the curve of L against m. For the fairly accurate approximate 
method of Timman, the maximum of m occurs slightly beyond separation, near which 
there is a fairly tight loop in the (L, m) curve. For the less accurate approximate 
method of Pohlhausen, the maximum of m occurs further beyond separation, and the 
loop in the (L, m) curve is a loose one. 


The near appearance of a cusp in the (L, m) curve for Timman’s method, and its 
actual appearance for the similarity solutions, suggests that it may be realistic to include 
it in the functions to be defined and, accordingly, to have the maximum value of m at 
separation. The functions H (m) and /(m) are chosen to satisfy these requirements 


FiGurE 8. The function /(m) according to several methods. 
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Ficure 9. The function H(m) according to several methods. 


and to join as smoothly as possible onto the values given by Thwaites (and amended 
by Curle and Skan) at separation. In the region downstream of separation, they are 
chosen to tend to the values for the similarity solutions as m—> 0. The function L (m) 
is chosen to be continuous with the amended Thwaites value, but with a cusp as 
described. 


The functions /(m), H (m), are given in Table I, while L (m), given approximately 
by L=0-45+6m upstream of separation, may be adequately approximated in the down- 
stream region by 

L=0:45+ 6m, 0:0900 => m > 0-:0625 
L=1-8m'/?+6m, 0:0625 >m>0 
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Ficure 10. The function L(m) according to several methods. 
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The Buckling of a Reinforced Circular Plate 
under Uniform Radial Thrust 


I. T. COOK and H. W. PARSONS 


(Department of Mathematics, The University, Hull) 


SUMMARY: An exact solution for the symmetrical buckling under uniform 
radial thrust is obtained for a thin circular plate having a particular type of 
thickness function for the cases in which the edge of the plate is either clamped 
or simply-supported. In both cases it is found that the critical thrust 
necessary to produce buckling can be increased from its value for the uniform 
circular plate of the same material and volume by concentrating material in 
the central region of the plate. For the clamped plate the increase is about 
18 per cent and for the simply-supported plate about 29 per cent. 


1. Introduction 


The initial symmetrical buckling of a uniform thin circular plate under uniform 
radial thrust applied round its edge has been discussed by a number of authors 
including Bryan“? and Timoshenko“. Meissner? discusses the problem of a 
uniform annulus and Olsson” the problem of an annulus with quadratically varying 
stiffness. 


In the present paper an exact solution is obtained for a thin circular plate of 
radius a having a simple thickness variation, buckling similarly under the same type 
of loading. Following Capey®’, the plate has a “step” thickness function (Fig. 1) 
so that it consists of an inner circular portion of constant thickness A, surrounded 
by an annular portion of constant thickness h,, the change-over radius from h, to 
h, being at. 

By varying h,/h, and t, the optimum distribution of a given volume of material 
in the form of a circular plate of given radius and having the thickness variation 
described can be determined for the loading and buckling pattern assumed. 


at 
FIGURE 1. 
Received March 1961. 
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NOTATION 
radius of plate 


h thickness of plate, — ~ eras 


p=h,/h, 
h_ average thickness of plate, =h, (1 — ft? + pt’) 
Poisson’s ratio (taken as 0-333) 

Eh’ 
12(1-«?) 


S parameter in the definition of T., 


D fiexural rigidity of plate, 


T applied radial thrust per unit length at circumference r=a 
12 (1 a? 


Six Maximum value of S for variable ¢ and p 


T. critical value of T for buckling, 


polar co-ordinates 


w=w(r), deflection of plate 


¢=dw/dr 

N,=B+A/r 

resultant stress components 
(A, B, C are constants) 
Nw=C/r’ 

do oo 
bending moment and shear stress 

ldo 

dr 

k?= B/D 

q=1-A/D 


D, A, B, C, k, q take suffix 0 or | according as O<r<at or at<r<a 
respectively. 


2. Equilibrium Equations for the Plate 

The stress resultants in the plane of the plate are determined by the boundary 
conditions at r=0, at and a and by the continuity of the radial displacement at r=at. 
It is found that N,=Ne=B,. Nw=0 for OSr<at and N,=B,+A,/r. 
Ne=B,-—A,/r’°, N.o=0 for at<=r=<a, where 


B, B A, T 
2p 
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BUCKLING OF A REINFORCED CIRCULAR PLATE 


Assuming that when the plate buckles it does so symmetrically, so that w= w(r), 
the equilibrium of a circular portion of the plate with centre O and radius r gives 


Substitution for Q and N, then gives the Bessel equation in ¢, namely 
1 do 


When 0 =r <at, q=1, and ¢=0 when r=0, so that solution of equation (3) gives 
(4,7) 0=r=<at, 
where E, is a constant. 
When ai = r<a, q=q, and (3) gives 
o=E,U (k,n)+F V (kyr), 
where E,, F, are constants and U (k,r), V(k,r) are independent solutions of 
equation (3). 
If g,=n’ > 0, then U (k,r) and V (k,r) can be replaced by the Bessel functions 


J, (kr) and J_, (k,r). It is to be noted that, in general, m is not an integer and the 
functions may be evaluated by the series expansion 


(k,r/2y" (k,r/2V (k,r/2y 


Jn (k,n) l’(n+1) (n+ 1) (n42) 


The gamma functions involved may be absorbed in the arbitrary constants 
E, and F,. 
If gq -m’ < 0, 
U (k,r) > [z, cos (m log k,r)+ 8, sin (m log k,r)} (k,r/2)* 


~ 


V (k,r) [y. cos (m log k,r)+ 4, sin (m log k,r)] (k,r/2)*, 


s=0 
mB,_, —mz2,_, —SB,_, 
where a,=1, 8,=0, «, , 8B, 
s(m? +57) (m? + s*) 


Y=0, %=1, s (m? +5?) s (m?* +s?) 


3. Determination of the Buckling Thrust 
E,, E, and F, must be chosen to satisfy the boundary conditions at r=dat, 
namely, ¢, M, and Q are continuous. Since, however, N, is continuous at r=at, 
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equation (2) shows that there are only two independent boundary conditions. The 
condition that ¢ is continuous at r= at gives 


EJ, (k,at)=E,U (k,at)+ FV (k,at) . (4) 


and the condition that M, is continuous at r=at gives 
DoE, | (koat)+ “J, (ka) D, (kat), + 


+ F, kV’ (k,at) + V (katy! |, (5) 
at 
primes (’) denoting differentiation. 


Using equation (4), equation (5) may be written 


(k,at)+ (p* —1) J, (K,at)| =E,k,U’ (k,at)+ F (k,at) (6) 
The boundary condition satisfied at r=a is that ¢=0 or M,=0, according as 
the plate is clamped or simply-supported. 


ie. E,U(k,a)+F,V (k,a)=0 (plate clamped) . (7) 
or 


E,|k,U’ (k,a)+ U (ka) +F, (k,a)+ V 0 (plate simply-supported). 
(8) 


Eliminating E,, E,, F, between equations (4), (6) and (7) or (8) gives the 
equation 


J, (k,at) U (k,at) V (k,at) 

k,atp’J,’ (k,at)+o (p* —1) J, (k,at) k,atU’ (k,at) kV’ (k,at) 
0 uk aU’ (k,a)+oU (k,a) uk,aV’ (k,a)+<V (k,a) 
0 (9) 


where »=0 or | according as the plate is clamped or simply-supported and k,, k, 
and q, are determined by equation (1), namely, 


p’k,” (1 -q,) 
2 2p-U-a)(p-1) ae D, 
The least positive root of equation (9) is 7.., the critical value of T for buckling. 
= Now the critical buckling thrust T for a uniform plate of radius a and thickness 
12 a?’ 
where A= 14-68 or 4-28 according as the plate is clamped or simply-supported. 


340 The Aeronautical Quarterly 


| 
4 
J 


BUCKLING OF A REINFORCED CIRCULAR PLATE 


TABLE I 
VALUES OF § FOR CLAMPED PLATE (FOR A UNIFORM PLATE, S= 14°68) 


t=003 01 02 03 O4 O45 O5 055 06 O7 O8 O99 097 


14°58 13:71 12°04 5°47 3-82 465 688 11:05 
0-4 9-14 
0°5 14:59 13°84 12°42 11°69 11:27 10°36 9°48 9-61 11-12 13-28 
06 12:07 12°60 13-89 
0:7 14°64 14:21 13:24 12°63 13-21 13:25 13-66 
0-9 14:66 14°56 14:25 13-94 14°35 14°57 1466 14°66 
14°68 14:77 15:00 15:39 15:37 15°16 15-00 14:43 14:45 
1:3 16°49 16°55 1645 16°16 15-71 13°53 
1-4 16:94 16°88 16°58 16:06 12°96 
14:70 14:91 15-49 16°94 17°15 17°16 16°91 16°38 12:99 12°39 
16 16°93 17°19 17:28 17:12 16°60 
7 16°79 17:07 17:24 17°16 
16°54 16°83 17°04 17-05 16°72 
19 16°23 16°70 16°79 
2:0 1469 14°83 15°13 15°85 16:07 16:27 16°40 16:32 14:85 12:03 10°09 11-32 
3-0 1466 14:45 13°73 11°58 10°80 10°62 11:22 11-11 8:26 7:87 
TABLE II 


VALUES OF § FOR SIMPLY-SUPPORTED PLATE (FOR A UNIFORM PLATE, S=4:28) 


0:2 0°6 0-7 0-75 0-8 0-85 0-9 0:97 


p=03 4:27 3°71 3-05 2°72 3:13 
0°5 4:27 3°78 314 2:96 4:08 
0°7 4:28 3:99 3°53 3:34 3°55 3:84 4°13 
0:8 3:98 
09 4:03 
4:29 4°34 4°47 4°58 4°58 4°53 4-42 
1:3 5-01 5:07 4:96 4°67 4:40 
4:29 4:39 4°70 5-17 5°35 5°36 5-29 4°88 4°47 
16 5:15 5°41 5°41 
5:07 5°40 5°50 5:49 5:33 4°52 
18 495 534 549 553 
1:9 4°79 5°54 
2-0 428 426 425 461 507 533 552 549 S21 457 
22 5:49 
3-60 401 4°89 5:33 5°33 

4:27 3-86 3:07 2°74 2:98 3°85 4°64 5°22 4-69 
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For comparison with T it is convenient to write 


SEh* 


12 (1-07) 


a and evaluate S. 


For given values of ¢ and p, equation (9) can be solved by trial and error. The 
value of t lies between 0 and 1, the extreme values of ¢ corresponding to a uniform 
plate for all values of p. The main features of the variation of S with t and p can 
be seen by considering values of p in the range 0< p< 3. Tables I and II show 
the results of the calculations for selected values of t and p. 


In Tables I and II, the heavy line separates the calculations for which the Bessel 
function expansions are used (above the line) from the calculations requiring the 
other expansions (below the line). From the tables, S,,..= 17:28 occurring for t=0-5 
and p=1-°6 (clamped plate) and Smax=5:54 occurring for t=0°8 and p=19 
(simply-supported plate), so that resistance to buckling by redistribution of the 
material forming a uniform plate can be increased by about 18 per cent or 29 per 
cent, according as the plate is clamped or simply-supported. 
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The Non-steady Irrotational Flow of an 
Inviscid, Incompressible Fluid, with Special 
Reference to Changes in Total Pressure 
through Flow Machines 


J. H. PRESTON 
(Department of Fluid Mechanics, University of Liverpool) 


SuMMARY: In this paper the way in which an inviscid incompressible fluid 

changes its total pressure is discussed. The basic equations are set out and it is 

shown that changes in total pressure are essentially associated with time 

variable flows. Simple examples are worked to illustrate this and finally 

application is made to single-stage flow machines, using a simple vortex 
representation. 


1. Introduction 


Dean“ has demonstrated that the frictionless flow through a flow machine, 
doing work or having work done on it by the fluid, must be unsteady. A little 
consideration shows that almost all practically useful flow machines exhibit this 
characteristic. As an example, consider a single-stage axial-flow machine such as a 
fan. The forces on the blades forming the fan do work on the fluid as it passes 
through the fan; this appears as a rise in total pressure of the fluid. The forces on 
the blades give rise to pressure differences on the two surfaces of the blades and in 
the spaces between the blades, and consequently there are associated variations in 
velocity. To an observer situated on the blade, the motion would appear to be 
steady; but to an observer at rest and near to the fan disc, the motion would appear 
to be unsteady and periodic, with variations in total and static pressure and in 
velocity and its direction. A simplified equivalent of the fan is the two-dimensional, 
infinite, moving row of aerofoils. If the whole system is given an equal, but oppo- 
site. velocity to that of the moving row, then the row is reduced to rest, becoming 
a fixed cascade of aerofoils, and the motion is steady. This is the basis of cascade 
testing and of calculating the forces and pressure distributions on the aerofoils. 
Since the aerofoils are at rest, they do no work and the total pressure is constant 
throughout the field of flow. If, at a fixed point near a moving cascade, information 
about the fluctuations in total pressure, velocity, and so on is required, then further 
investigation is necessary. So far as the author is aware this has not been previously 
attempted and in the past has been deliberately avoided by use of the concept of an 
actuator disc. 
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In this paper the basic equations for the unsteady irrotational motion of an 
inviscid incompressible fluid are set out. Simple, illustrative applications of these 
are made (a) to the problem of the flow along a long pipe connected to a large 
reservoir when a valve is suddenly opened and (5) to the fluctuations in total 
pressure at a point when a vortex or a doublet are in uniform motion past this point. 
The calculations are then extended to an infinite moving row of bound vortices and 
the variation of total pressure with time at a fixed point is obtained. This gives a 
simplified picture of the total pressure changes through an axial flow fan. The 
problem of a single-stage radial-flow machine is similarly approximated by a 
rotating ring of equally spaced bound vortices and the solution for the infinite 
straight row of vortices can be extended to this problem by a simple conformal 
transformation. Calculations have also been made for a moving row of doublets to 
indicate the effect of aerofoil thickness. The actuator disc concept is examined on 
the basis of the total pressure fluctuations. Also, interaction of the fixed and moving 
stages in a compressor is briefly discussed. A short discussion is also given of the 
significance of the results in relation to the design and calibration of instruments for 
measuring total pressure, static pressure, velocity and flow direction in tests on 
flow machines. 


NOTATION 
q_ resultant velocity 
s distance measured along a streamline 
D/Dt_ differentiation following a fluid element 
t time 
p Static pressure 
P total pressure 
p density (constant) 
g acceleration due to gravity 
z height above some datum 
velocity potential 
P, datum value of P in time-variable flows 
H_ height of fluid in reservoir 
length of pipe 
U_ velocity of flow in pipe (in Section 3.1 only) 
Pp, atmospheric pressure 
P, total pressure at pipe exit 
W potential function 
z=x-+iy, complex co-ordinate 


I’ strength of bound vortex (circulation) 
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u,v velocity components 
strength of doublet 
U_ velocity in x-direction of moving vortex or doublet 
V_ velocity in y-direction at infinity 
L, x-component of force on bound vortex 
L, y-component of force on bound vortex 
h_ pitch of vortices or doublets in a row 


n number of vortices 


P, total pressure at y= —°o 

P, total pressure at y= +00 

U, x-velocity component at y= 
U, x-velocity component at y= +o 
AP=P,—P, 


(2 angular velocity of rotating ring of vortices 
a radius of ring 
L, force component in circumferential (@) direction 
m source strength 
v_, radial component of velocity at radius a 
C= E+in 
- 4+ 16>complex co-ordinate in plane of ring of vortices 
2X length on real axis of z-plane representing perimeter of circle radius a 
in (-plane. 


bar, denoting time average over a cycle 


2. Basic Theory 


Consider the motion of a fluid element along a streamline. Denote differentia- 
tions following the fluid element by D/ Dr; then, if g is the velocity and s is measured 
along the streamline, 


(1) 

Dt ot as° 
Denote the total pressure by P, then 

(2) 


where p is the static pressure, p is the density and z is the height of the element 
above some datum. 
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The equation of motion along the streamline is 


1 oP 


= 3) 
ot p os 


This will form the basis of the subsequent calculations. The change of total pressure 
of a fluid element is 


DP oP oP 
= +q 
Dt ot os 


The second bracket term is zero, from (3). 
DP 
Hence (p + pgz) 
Op 


Hence the total pressure P of a fluid element only varies with the time if the static 
pressure at a fixed point varies with the time. Since the purpose of a flow machine 
is to give or abstract total pressure to or from the fluid, the motion relative to a fixed 
point must be unsteady. 


The result in equation (4) is contained in a similar result for stagnation 
enthalpy derived by Dean‘ and it is apparently given in a book by Eck’. 


Returning to equation (3), since the motion is irrotational a velocity potential 
® exists and 


Co 

(5) 
whence, integrating (3) along a streamline, 
Co P 

+ =f (t). ; 6 

ct Pp (6) 


This equation is given in Lamb’s Hydrodynamics (paragraph 20). 


3. Applications of the Theory 


3.1. LONG UNIFORM PIPE CONNECTED TO A LARGE RESERVOIR. SUDDEN OPENING 
OF VALVE AT END OF PIPE 


This is a textbook example but it is given here to show the part which 
acceleration plays in determining the total pressure. 
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PRESSURE CHANGES IN FLOW MACHINES 


Po 
H 
= 
U 
E 
(a) 
—— 
V(2gH) (b) 
0-6 — + — 
c) — — 
02 


1-O 2-0 3:0 


tv(@gH) 


Ficure 1. Flow from a reservoir. 
(a) Notation. 
(b) Velocity against time (equation (10)). 
(c) Loss of total pressure against time (equation (11)) 


The level of the liquid in the reservoir above the axis of the pipe is denoted by 
H and is assumed to remain constant (see Fig. 1(a)). 


Equation (6) therefore becomes 


> + P = constant = P, (7) 


which holds throughout the fluid, since the motion is irrotational and the fluid is at 
rest in the reservoir. Denote the length of pipe by L, the velocity in the pipe by U 
and atmospheric pressure by p,, which is applied at the reservoir surface and at the 
pipe exit after the valve opens. Contributions to ¢ in the neighbourhood of the 
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pipe entry are neglected and the value of ¢ there is taken to be zero. Applying (7) 
to stations at the entry and exit of the pipe, we have, if p, is the pressure just inside 
the entry, 


L 


p+ U dr, 


0 


which is the equation of motion for the fluid in the pipe. Between the pipe entry 
and the reservoir surface, , and hence 0¢@/0t, is negligible and so 


p,+4pU* =p, + pgH, (9) 
whence (8) becomes 
pgH —4pU? = pL 
or dt=2L 


The condition t=0, U=0 gives C=0 


U t/ (22H) 
( 
and /QeH) tanh . . (10) 


This is shown plotted in Fig. 1(b), where it is seen that U rises from U =0, t=0 to 
U= /(2gH), t=0o. 


Denote the total pressure at the reservoir surface by P, and at the pipe exit by 
P, then 


U? 
P, — Po= p> 
P,—Pr 


t/ (22H) 
pgH 2¢H 


= sech? 


This is shown plotted on Fig. l(c). It will be noted that, at t=0, U=O, 
P,—Px=pgH, so that the available head H is used in accelerating the fluid in the 


pipe from rest. At t=00, U= /(2gH) and (P,—P,)=0, since the motion has now 
become steady. 
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If the level of the reservoir falls perceptibly or is controlled in some way as a 
function of the time, then the “constant” of equation (7) must be replaced by the 
appropriate function of the time. 


One further point is worthy of note: 


Hence (p, —~p,)dt=pLU(t). . 


0 


Hence the “pressure-time’’ integral will yield the velocity or mass flow. Practical 
use is made of this in the Gibson method of determining the mass flow through the 
penstock of a hydraulic turbine installation. 


4. Vortex Moving Past a Fixed Point 


il 
W= o+i¥ = logz. 


gives the potential for a stationary vortex at the origin. 


(14) 
6=tan-' (y/x) 


Suppose the vortex is in steady motion with velocity —U parallel to the x-axis. 
Then the co-ordinates of a fixed point relative to the vortex are x= Ut and y, 


ct Cx 
dw 
dz 
whence, with x= Ut, 
y 
“=~ (Un +y?’ 
iy 


2x (Ut)? +y’ 
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(a) 


-0-2 


Figure 2. Uniform motion of vortex and doublet pass a fixed point. Change of total 
pressure with time. 

(a) Vortex (equation (19)). 

(b) Doublet (equation (23)). 


Equation (6) may be written 


0 + 


p +P=constant=P,,_ . ; 8) 


since the motion is irrotational and at rest at infinity. Hence, using (15) and (16), 


4 y 
Fs P 4 Qn (Ut)? 
Ty= 9 


This is shown plotted in Fig. 2(a). The velocities are given by (16) and (17), so 
that the direction of the velocity vector and the static pressure can be found as 
functions of the time. 


5. Doublet Moving Past a Fixed Point 


For a stationary doublet of strength 27v with its axis parallel to the x-axis 


v 
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2xy 

22 

( x? ( ) 


For a doublet of sirength 2<v moving with a velocity —U parallel to the x-axis 
(which is the same as a circular cylinder of radius /(v/U) moving with this 
velocity), 
x=Ut 
and, from equations (15), (18) and (21), 
= pUv (Ut/y)-1 
at y? {(Ut/y) +1}? 
(P-—P,) y’ 


This is shown plotted on Fig. 2(b). 


6. Moving Row of Bound Vortices 


Consider first the flow arising from a fixed row of vortices (Fig. 3(a)). Let the 
vortices lie along the x-axis with a spacing A and a circulation |" in the anticlock- 
wise direction. 

The velocity potential W (Lamb paragraph 156, or by use of the transforma- 
tion (37)) for the vortices only is 


and the velocity components are 


—T sinh (2zy/h) (25) 
2h cosh (2zv/h)—cos (27x/h) 
4 sin (27x / 
and 


2h cosh (27y/h)—cos (27x/h)° 


Consider now conditions at a fixed point when the row is in motion. Then, 
relative to instantaneous axes moving with the row, the fixed point has co-ordinates 


x=Ut and y, 


which, inserted in equations (25) and (26), give u and wv at the point. 
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" Figure 3. Vortex row or vortex ring in uniform motion. Total pressure changes. 
(a) Row of vortices in uniform motion. 

(b) Total pressure against position. 

(c) Rotating ring of vortices. 
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Ficure 3(d). Total pressure against time 
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To calculate the changes in total pressure, we use equation (6) in the form 
p +P=constant=P,,_ . . 127) 


since the motion is irrotational and steady at y= +00. 


With x=Ut and using equation (25), 


Co 
P,—P=pe pUu 


_ —pUr’ sinh (27y/h) (28) 
2h cosh (2y/h)—cos (27Ut/h) 
For y= —© and y= + © respectively the total pressures are 
_p —p — eur 
| 
so that the gain in total pressure is 


From equation (28), the variation in total pressure is given in Fig. 3(b), where 

2h(P—P,)/pUV has been plotted as a function of 2~y/h for a range of values of 

Ut/h. It may be noted that the fluctuations in total pressure have 

practically vanished at y/h=1-0 and amount to 10 per cent of the total pressure 

change through the row at y/h=4. A cross-plotting in Fig. 3(d) shows the variation 
4 through the cycle 0 < Ut/h < 1-0 at constant values of 27y/h. 

The average value of (P—P,) over a cycle (0 << Ut/h < 1:0) can be computed. 


From equation (27) 


where the bar denotes the mean value over a time T=h/U. 


(P- P.) — (ur h=9 
From equation (24) it is readily found that 


ih 
tan~' [(cot (tUt/h) tanh (zy/h)). ; 
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Hence, for y positive, 


(P—-P,) UU 
p 2h | 


and for y negative, (33) 


(P—P,)_ | 
p 2h 


So that the average value of (P—P.,) is independent of y on either side of the vortex 
row, but it changes discontinuously through the row by |'U/h, thus substantiating 
the actuator disc hypothesis. 


%7. Rotating Ring of Bound Vortices 


n bound vortices of strength |’ are disposed at equal angular intervals at a 
radius a and rotate with angular velocity {2 (see Fig. 3(c)). There is a radial flow 
through the rotor which can be represented by a source of strength m at the origin. 


The rate of doing work by the rotor is 
x opposing torque={2x —nL,a, 
where L, is the force on the vortex in the direction of @ positive, 


m 


Lo p —pv, T, (34) 


where 7, is the radial component of velocity at r=a. Hence the rate of doing work 
equals 


(Pm). 
This must equal the rate at which the fluid gains total pressure (AP), i.e. APm. 


Hence AP a (35) 
From equation (31), 
P-P,=-p (36) 
ot 


so that it is necessary to determine ¢ for a stationary ring of vortices. 
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This is done by transforming conformally the ring of vortices into an infinite row 
of vortices. 


where (=£+in is the plane of the ring and z=x+iy is the plane of the row. 
Putting (=ae*", ‘ ; 
so that r= /(£?+7n*)=ae, 


then p=—ty/X 
6=n7x/X 


(39) 


A strip of width 2X, stretching from y=-—0°O in the z-plane, maps 
into the whole of the (-plane and, in particular, the portion of the x-axis 
0<x<2X corresponds to the full circle radius a and 0<.4@<2r. Also, the 
vortices which are disposed at equal intervals of 6=2z7/n on the circle r=a (u=0) 
in the (-plane correspond to equally spaced vortices (pitch A) on the x-axis in the 


z-plane. 
2x th 
Using equation (39), 
so that ‘ ‘ . (40) 


Since successive strips of width 2X in the z-plane map into similar successive 
sheets in the (-plane, we have therefore transformed a ring of vortices in the (-plane 
into an infinite row of vortices on the x-axis in the z-plane for which, from (24), 


il TZ 
W=o+i¥=— logsin —~ 
It therefore follows that the values of ¢ are the same at corresponding points in the 
two planes. 


When the ring of vortices rotates with angular velocity 2, then, for a fixed 
point on the €-axis, the co-ordinates relative to the rotating system are « and 
6= —(1; but, from (39), 


p= —ty/X 
O=7x/X 
U 
and, with x=Ut, Q=-« a (41) 
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Using relations (39), (40), and (41), the value of é@/ét (which is the same at 
corresponding points) can be found from equation (25) for the moving row of 
vortices and thus 


poe 
P-P,=- 
ot 
sinh nu 
4x coshnu—cos nQt’ 
(P—P.,) sinh nu 
or (42) 


prQ on cosh nu—cos nQt 


which has the same form as equation (25) and so the Figs. 3(b) and 3(d) still apply 
with appropriate change of variables. The change in total pressure from (“= 00, 
to (u= + 00, r=00) is 


which is identical with equation (35). 


8. Moving Row of Doublets 


To complete the picture the calculations are now made for a moving row of 
doublets which give a first approximation to the effects of thickness. Since there 
is no force on the doublets they do no work and so there is no net gain of total 
pressure across the row. 


The velocity potential W for a row of doublets of strength v and pitch A situated 
on the x-axis with their axes parallel to the x-axis can be obtained by differentiating 
equation (24) or using the transformation (37) to transform a doublet and its image 
in the half ¢-plane into an infinite row in the z-plane. 


sin (27x/h) 


h [cosh (2zy/h)—cos (2x/h)] 


This yields 


‘=\? [1 —cos (27x/h) cosh (2zy/h)] 
u eV h [cosh (2=y/h)—cos (2=x/hA)P (45) 
=—2v \* sin (2=x/h) sinh (2=y/h) 
2 [cosh (27y/h)—cos (46) 
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The variation in total pressure for a moving row is found as before by putting 
x=Ut. Then 


Op 
P-P,= ot —pUu 
which, using equation (45), yields 
P—P, 1 Uh?  [l—cos Ut/h) cosh (2ty/h)] (47) 


pU? [cosh(2zy/h)—cos (2% Ut/h)} 


This is shown plotted in Fig. 4. 


Over a cycle 0< Ut/h <1-0, the average value of 0¢/¢t is, for the period 
T=h/U, 


U 
=[(9)utn=, —@vtn=ol h 
whence, using equation (44), 


and the average value of the total pressure is P,. 


9. Concluding Remarks 


The first example (Section 3.1) shows that, in the accelerating flow from a 
reservoir, the total pressure is not constant. The assumption of an incompressible 
fluid means that, on the opening of the valve, the pressure adjusts itself instan- 
taneously along the pipe for a constant value of (p,+pgH) to a value varying 
linearly from (p, + pgH) at the inlet to p, at exit. In a real fluid the adjustment will 
take place over a finite time with a pressure wave travelling up the pipe at the speed 
of sound and with possible reflections at each end. The use of a pressure-time 
integral to obtain the mass flow may be noted. The time delay which occurs before 
steady mass flow conditions are reached after the sudden opening of a valve may 
have an application to the fuel control of aircraft engines. The set-up might 
possibly be used to check the response of instruments designed to measure total 
or static pressure. 


The motions of a vortex and of a doublet (given in Sections 4 and 5) relative to 
a fixed point set up time variations in total pressure, static pressure and velocity at 
a fixed point. In combination they could be used to represent the effects of the 
passage of a lifting aerofoil of finite thickness as experienced by a distant point. 
Again, it is possible that the solutions might be of use in an experimental check of 
the response of measuring instruments. 


Sections 6 and 7 are concerned with effects experienced at a fixed point caused 
by a moving row of bound vortices and a rotating ring of bound vortices. They 
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(a)TOTAL PRESSURE AGAINST POSITION 


Ficure 4. Row of doublets in uniform motion. 
(a) Total pressure against position. 
(b) Total pressure against time. 


give a first approximation to the effects of the blades of single-stage, axial flow or 
radial flow machines. In both case the time variations in total pressure at a fixed 
point are entirely similar. The variation is periodic and is anti-symmetrical for 
equidistant points on either side of the moving row. The variations of total pressure 
close to the row are large, but fall away to steady values at large distances. Ata 
distance about equal to the pitch of the vortices the fluctuation in total pressure is 
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less than | per cent of the total pressure rise through the row, and at half the pitch 
it is about 10 per cent. Since, in practice, the distance between the rotors and 
stators of an axial flow compressor is about half the pitch, instruments for measur- 
ing total pressure, and so on, will be situated at about a quarter of the pitch from 
the moving row. They will thus be subject, on the basis of the present theory, to 
considerable fluctuations, which may cause errors in the mean readings. It may be 
noted that the time average of the total pressure is constant on either side of the 
vortex row, but it changes discontinuously through the row, thus agreeing with the 
actuator disc concept. 


If we consider a stator in front of or behind a rotor, then it will be subject to 
fluctuations in magnitude and direction of the velocity. Consequently the forces on 
the stator blades will change with time and free vortices will be shed downstream 
and thus the blades of the rotor will find themselves passing through a fluctuating 
flow field. A complex interaction will be set up between the stator and rotor, giving 
rise to fluctuating loads on both the stator and rotor blades, which may be of 
significance in relation to fatigue and flutter. 


The present calculations could easily be extended to cover continuous 
distributions of vorticity and doublet strength along the aerofoil chords and also if, 
for aerofoils in a fixed cascade, the velocity field is calculable, then the time 
variations in total pressure, and so on, at a fixed point in relation to a moving 
cascade, are also calculable. 


10. Conclusions 


A basic theory exists for the calculation of the total pressure changes at a fixed 
point in the non-steady, irrotational motion of an inviscid, incompressible fluid. 


This has been applied to a number of problems, in particular to the calculation 
of total pressure changes with time at a fixed point in relation to the flow through 
a simplified single-stage flow machine. The periodic variation of total pressure 
with time varies in magnitude and form with the position of the fixed point and it is 
found that fluctuations in total pressure are less than | per cent of the total pressure 
rise through the stage at distances greater than the pitch of the vortices. 
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The Deformational Behaviour of Low Aspect 
Ratio Multi-Web Wings* 


ParT I: Experimental Data 


R. H. GALLAGHER? and I. RATTINGER?t 


SUMMARY: Results of a study of the accuracy attainable from various 
approaches to low aspect ratio wing deformational analyses are described. 
Seven model multi-web wings, representing a restricted range of sweep angles, 
aspect ratios and other geometric factors, were tested for deflection influence 
coefficients. Methods of matrix structural analysis and other techniques, 
involving numerous idealisations, were then applied in the prediction of the 
test data. This paper presents the experimental data. Subsequent papers will 
describe the analytical methods, the level of test-theory correspondence, and 
the merits and shortcomings of the t.eories studied. 


1. Introduction 


In the period following the Second World War aerodynamic efficiency 
considerations dictated a change from wings of essentially rectangular and high 
aspect ratio plan form to swept, low aspect ratio plan forms. Continued progress 
in design resulted in still other types of plan forms (e.g. the delta wing) as well as a 
trend towards wings of high solidity. Analysts responsible for the prediction of the 
elastic characteristics of proposed supersonic wing designs recognised the need to 
use more sophisticated concepts of structural behaviour than those encompassed by 
elementary beam and torsion theory. Both indeterminate structural analysis pro- 
cedures and approaches based on classical plate theory found favour. Whatever 
approach was taken, however, the analyses involved the solution of higher order 
simultaneous equations than could be solved by hand computation. Fortunately, 
a concurrent development of automatic computational devices, both analogue and 
digital, occurred during the post-war period. The devices made feasible the solution 
of problems phrased in terms of many unknowns and, in addition, it has become 
economical to formulate automatically the equations to be solved. A major 


*The contents of this paper formed a part of a presentation to a meeting of the Structures and 
Materials Panel, Advisory Group for Aeronautical Research and Development, N.A.T.O., in 
Athens, 28th September 1960. 

+Aerospace Department, Bell Aerosystems Company. 

tAerospace Corporation. 
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remaining problem lies in the selection, from the many proposed theoretical 
approaches, of that one which will accurately and efficiently predict the elastic 
characteristics of a given wing. The efforts of the study described in this paper 
were specifically directed towards answering questions of accuracy. 


Two primary factors upon which accuracy depends are the analytical 
representation of the structure and the number of analytical degrees of freedom 
employed. Although it should be apparent that greater accuracy will be obtained 
when more realistic analytical representations of the structure are used, this factor 
has often been obscured. Questions concerning the relationship between attainable 
accuracy and the number of analytical degrees of freedom have not been considered 
to a significant extent in the literature. To evaluate the relationship between these 
factors, deflection influence coefficient tests were performed on seven realistic, yet 
simplified, low aspect ratio wing models and the results were compared with 
analytically predicted values. Each of the wings was of full span but was tested 
only for symmetrical loadings. The models were chosen so as to introduce, in a 
controlled manner, those problem variables that could be studied within the limited 
number decided upon for the programme. 


2. Experimental Programme 


Deflection influence coefficient tests were performed on seven model wings 
designed specifically for the purpose of delineating the effects of sweep angle, aspect 
ratio, relative skin thickness, and internal member configuration. Other significant 
factors such as plan-form taper and variations in the type of internal structure and 
support conditions were represented in the three additional wings not discussed in 


CENTRE 2 
UNE AXIS 
21N. WIDE SUPPORT ~ LIN. Zw 
Q 
| = 
y 
| 
i 


— 
PLAN FORM ad 
FACE OF SUPPORT .| FACE OF SUPPORT 


2 IN. WIDE MACHINED SECTION 8-8 
FILLER 8LOCKS 


Ficure 1. Model 1. Geometry and Getails. 
Model material: 6061-T6 aluminium. All internal members formed of 0-040 in. sheet. 
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FiGuRE 2. Geometry and node point locations. Model wings. 


All models are clamped along centre support. O denotes a node point. 


this paper. Data relating to these three wings, as well as more detailed information 


on the tests to be described, are given in Ref. 1. 

Parametric variations of the foregoing factors essentially have their basis in 
the rectangular plan form, aspect ratio 4-0 wing (Model 1) shown in Fig. 1. This 
and all other models were fabricated from 6061-T6 aluminium alloy sheet. All 
interna' members were formed channels, 0-040 in. thick, while the skins were single 
sheets, 0:063 in. thick. The spars were 30 in. long for each semi-span and 
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continuous from tip to tip, being joined to the rib segments with clip angles. The 
skins were joined to the internal members with } in. diameter rivets, spaced } in. 
on centre. 


During test the model was clamped in a 2 in. wide support along the X-X axis 
(see Fig. 1). Severe clamping effects were precluded by means of a set of solid 
internal blocks in the support area. Load application points were defined by the 
spar-rib intersections, numbered as shown in Fig. 2. The remaining six model wings 
were of the same basic construction, with the exceptions noted in the following 
descriptions (see Fig. 2). 


Models 2 and 3 were swept 45°; the root-to-tip spar lengths of Model 1 were 
retained. Model 3 had a skin thickness of 0-125 in., and was the only model in the 
series to differ in this regard. Model 4 was unswept and of half the span of 
Model 1. Model 5 was swept 30°, and semi-span spar lengths of 30 in. were again 
retained. Models 6 and 7 were delta wings, the former representing the 45° swept 
wing (Model 2) with structure added to form an unswept trailing edge. Model 7 
typified the rectangular plan form wing (Model 1) with an added 45° swept leading 
edge area. Note however, that the tip chords and spans of these models differ in 
length from those of Model 1. 


The manner in which an ordered variation in the chosen parameters was 
effected is now discussed (see Fig. 2). Models 1 and 4 attempt to highlight the 
effects of aspect ratio, chordwise bending, proximity of load to the support and shear 
lag. Sweep angle effects were delineated by Models 1, 2, and 5, where there is 
consistency of skin thickness, depth of section and internal member configuration. 
Models 2 and 3 assess the effects of skin thickness, while internal structural 
configuration effects were defined in Models 6 and 7. 


Typical test arrangements are pictured in Fig. 3, which shows Model 6 
positioned in the support frame. To ensure symmetric behaviour about the X-X 
axis, loads were applied to corresponding points on both semi-spans through rods 
with hemispherical ends. The loading apparatus was supported independently of 
the model support fixture. Load variations were accomplished by a turnbuckle 
arrangement. Load magnitudes were measured by a 1000 Ib. Baldwin load cell 
placed in series with the rods and deflections were measured directly from dial 
gauges that could be read to an accuracy of 0-001 in. at points near the support and 
0-01 in. at outboard points. 


Before testing, a number of preliminary operations were conducted in order to 
minimise experimental error. All points were cyclically loaded to remove initial 
slippage, to define the extent of hysteresis and to aid in determining the load range 
and intervals to be used during test. It should be noted that hysteresis effects proved 
to be insignificant for all models. 


The actual test sequence started by recording all zero-load dial gauge readings 
Then, the load was increased to one interval above the maximum value desired and, 
as the load was removed in decrements, deflection measurements were made at each 
level. Complete unloading was comprised of six equal decrements. Upon com- 
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FiGurE 3. Typical model wing testing arrangement. 


pleting this sequence the initial and final zero-load dial gauge readings were 
compared. Re-testing was undertaken if variations in the deflection differences 
between successive load decrements or deflectional dis-symmetry about the X-X 
axis were greater than 10 per cent, or if there were indications of permanent set as 
given by zero-load gauge reading comparisons. 


Once a satisfactory set of readings had been acquired, deflection influence 
coefficients were obtained from the slope of a least-square straight line fitted through 
the load deflection points. Slopes were calculated on the basis of not only all load 
intervals but also all intervals except the lowest, since any existing hysteresis effects 
were most pronounced near zero load. A careful study of the measurements was 
made to determine the most representative value and this was the last step in 
reducing the test data. 


3. Results 


Due to the stringency of the criteria for the suitability of recorded data the 
cross-coefficients in each complete set of displacement influence coefficients ([4]) 
were closely in accord with the reciprocity condition. This condition requires that 
8,,=5,, where 5, symbolises the displacement at point i due to a load at point j and 


8, is the displacement at j due to a load at i. (All loads and displacements discussed 


November 1961 365 


» ‘ 


JAPOW- Uo 


SI vl I ol 6 8 L 9 € I 
6910-0 | L$70-0 1800-0 OF00-0 8000-0 6070-0 810-0 7770-0 7800-0 LLE0-0 670-0 $170-0 
8110-0 | LEI0-0 $800-0 600-0 $7Z00-0 8hZ0-0 8170-0 SLI0-0 610-0 £600:0 Sst0-0 61€0-0 CLZ0-0 0610-0 vl 
Z 200-0 1800-0 $010-0 | 8010-0 I[200-0 9600-0 0910-0 IS10-0 StI0-0 8010-0 St70-0 OF70-0 6070-0 6810-0 ZZL10-0 tl 
700-0 900-0 7800-0 9600-0 £800-0 8800-0 9600-0 9710-0 6210-0 £€Z10-0 O£€10-0 S€I0-0 9510-0 Ti 
8700-0 800-0 1200-0 710-0 ¥L10-0|LL10-0 8100-0 6£00-0 $$00-0 0010-0 LrI0-0 7200-0 9900-0 0600-0 8710-0 
9870-0 $070-0 OL10-0 SOI0-0 | LOOI-0 6180-0 6S£0:0 6070-:0 660-0 1840-0 Ol 
6070-0 9610-0 8210-0 8510-0 6£10-0 £850-0 0850-0 | L790-0 79P0-0 LEEO-O SIZO-0 660-0 St80-0 9190-0 | 6 
™ $910-0 ZLI0-0 O810-0 8910-0 7Z6h0-0 SIS0-0 | 87£0-0 6£80-0 0940-0 190-0 7Z8S0-0 Z0S0-0 8 
Q Of10-0 ISI0-0 £€L10-0 8210-0 8070-0 9760-0 8960-0 ZISO-0 | ¥SE0-0 8870-0 78S0-0 8rS0-0 LZ7S0-0 S0SO-0 | L 
9010-0 9E10-0 €L10-0 9170-0 OSZ0-0 6LE0-0 7Z0SO-0 1650-0 1890-0 | 6££0-0 ZELO-O ESEO-0 HLEO-0 | 9 
810-0 7670-0 7970-0 ¥EZ0-0 #810-0 FIOI-0 LP80-0 9920-0 9020-0 081-0 | SILI-O 6IET-0 8960-0 | 
870-0 9170-0 970-0 9170-0 160-0 $680-0 $80-0 080-0 991-0 SP9l-0 | 761-0 LOZI-O0 LI60-0 
LSZ0-0 £780-0 PS80-0 9£80-0 TZE80-0 99ST-O CLSI-O O9ET-O LITI-O 680-0 
© ZTIZ0-0 8€70-0 LS%70-0 6470-0 0870-0 LrL0-0 7080-0 9680-0 0880-0 160-0 PSST-O IE9T-O | 6660-0 £980-0 
SL10-0 9170-0 ZSZ70-0 $870-0 ZI£O-0 0990-0 0S40-0 780-0 7060-0 1860-0 9671-0 SIPI-0 OIST-O HS80-0 I 
(sjapou 
vl tl Ol 6 8 L 9 € I 410q) 
sjuiod 
x | uO siulog uonvayddy 
pooT 
= 
(Speo] 0} ONp ‘soyou! UI UMOYS 2 
Z GNV | STAGON JONINTANI 
513008 


t 


LOW ASPECT RATIO WING DEFORMATION 


7600- 
OS00- 
8700- 
1100: 
£810- 
6710: 
$800: 
£S00- 
9670: 
L970: 
917Z0- 
S9l0- 
710: 


SI 


0 


vi 


0 | 1800-0 9600-0 


8800 


8100: 
OS10- 
8010- 


0800 - 


LS00- 
81Z0- 
7610- 
9910: 
6010- 


vl 


-0 | 9500-0 


0 900-0 
0 sv00-0 
0 O£00-0 
0 1010-0 
0 600-0 
0 0600-0 
0 1800-0 
0 4900-0 
0 L¥I0-0 
0 8st10-0 
0 8710-0 
0 
0 010-0 


tl 


tl 


9000: 
€100- 
€Z00: 


6900-0 |.0900-0_ 
£500-0 6800: 
9500-0 $Z00- 
0900-0 
$900:0 
7200-0 8500: 
9100-0 $800- 
0800-0 1£00- 
6800-0 ¢P00- 
0600:0 7ZS00- 
€600:0 $900: 


0 
0 
0 


0 


ol 


7000-0 7600: 
9000-0 £900: 
9100-0 Z£00- 
6100: 
| 7600-0 8000: 
1LS0-0 | Z810- 
SOr0-0 
8820-0 
9020-0 1610: 
0Z10-0 9710: 
1090- 
1890-0 1€S0- 
rrs0-0 
£9€0- 


$600:0 8200-0 ZLIt0-0 


Ol 6 


0 
0 
0 
0 
0 


0 
‘0 £€10-0 600-0 
-0 $€Z0-0 | Zz10-0 


/apopy uo 


6 8 


0900-0 
600-0 100-0 
$£00-0 9600-0 
9700-:0 OF00-0 
9100-0 
££10-0 $800-0 


L 9 


000-0 
8Z00°0 0700-0 
6£00:0 ££00-0 
9900-0 8900-0 
¥900°0 010-0 
O£00-0 
8900-0 $00-0 
8600-0 2800-0 


0 0610-0 6810-0 |sst0-0 


0 
0 
0 
0 
0 
0 


9S10-0 
OSt0-0 
800-0 
9L£0-0 $870-0 
970-0 
7L7O-0 1870-0 


8 L 


uo SJulod 


1810-0 | Z0z0-0 
1610-0 SLbl-0 
1610-0 OSII-0 
2020-0 0060-0 
8020-0 ZIL0-0 


v 


0010-0 0400-0 
9900-0 7$00-0 
9700-0 0S00-0 9500-0 
0£00-0 OF00-0 ¥S00-0 
8700-0 600-0 


O810-0 8 


0810-0 9210-0 6£€10-0 
9FIO-0 7910-0 
7600-0 6110-0 OSI0-0 
$900-:0 4600-0 
800-0 6720-0 
8860-0 | 7S7Z0-0 
7180-0 9740-0 | 0170-0 
8850-0 
1€Z70-0 L7S0-0 68b0-0 £960-0 


t 


Il 


FININTANI NOLLO 


VL 


41490 LSAL 


7600-0 
9110-0 
LS10-0 
610-0 
8810-0 
8910-0 
8070-0 


600-0 
19¢0- 


N 


I 
0100-0 
Z£00-0 
0S00-0 
9800-0 
O110-0 
1900-0 | Ol 
£600-0 6 
| 
1610-0 | L 
9670-0 9 
9110-0 | 
OL10-0 | 
| € 
L7€0-0 | 
LOvO-0 | 
yjoq) 
| syuiod 
poo] 
7300W 
T300W 


| 
| 
| 
| | | 
| 
| 
| | | 
| 
— 
| | 
| | 


L /2POW uo S]U10d 


8Z10- 0| 6810- 0 9600-0 0F00-0 8100-0 0100-0 1000-0 £000-0 7910-0 OL10-0 2900-0 6£00-0 8100-0 7000-0 L910-0 SZ10-0 9800-0 9S$00-0 T£00-0 81 
¥800-0 ££10-0) 8010-0 8F00-0 9700-0 F100-0 2000-0 7000-0 €710-0 9200-0 6700-0 0£00-0 7100-0 LZ710-0 9600-0 8900-0 800-0 Ll 
1$00-0 7800-0 €810-0) #900-0 S€00-0 1700-0 7100-0 9000-0 400-0 7200-0 100-0 8500-0 7700-0 9700-0 8600-0 F600-0 8800-0 LL00-0 0900-0 91 
¥S00-0 8810-0| €200-0 0Z00-0 1100-0 OF00-0 6700-0 7900-0 7100-0 0900-0 $900-0 0100-0 9200-0 1800-0 7200-0 SI 
0700-0 6700-0 ZS00-0 4600-0 €£10-0| ¥L00-0 LE00-0 0700-0 700-0 I1£00-0 £P00-0 0900-0 7200-0 6S00-0 1S00-0 0900-0 £200-0 0800-0 vl 
7100-0 1700-0 ZP00-0 $900-0 1600-0 7800-0 OF00-0 9000-0 L100-0 700-0 LE00-0 7200-0 9100-0 9700-0 9£00-0 0S00-0 9900-0 tl 
¥000-0 0100-0 £700-0 OF00-0 £S00-0 8010-0 I9I0- -0/€110-0 0 7000-0 9000-0 0100-0 7Z00-0 8£00-0 0800-0 1000-0 8000-0 0700-0 Z£00-0 O0S00-0 Zl 
6910-0 8710-0 0010-0 9200-0 1$00-0 I1£00-0 0100-0 68€0-0 7970-0 $910-0 €O10-0 7900-0 $Z00-0 FHTO-0 7910-0 £010-0 II 
ZZ10-6 9FIO-0 OZIO-0 8400-0 SS00-0 8700-0 870-0 IPE0- 0| $9z0-0 €Z10-0 ZI10-0 6200-0 9600-0 O8f0-0 6210-0 7E10-0 Ol 
£800-0 8600-0 8S10-0 OZ10-0 0800-0 €F00-0 9170-0 9170-0 -0|9610-¢ 0 P£10-0 8600-0 7900-0 8SZ0-0 0970-0 L¥Z0-0 9070-0 0910-0 6 
7900-0 ¥L00-0 6010-0 bS10-0 FIIO-0 0200-0 H910-0 Z0Z0-0 $970-0 €810-0 8£10-0 7600-0 ZLI0-0 9810-0 £070-0 8IZ0-0 S610-0 8 
0v00-0 900-0 8800-0 OFIO-0 8FI0-0 9010-0 8110-0 8S10-0 F070-0 9970-0 87Z£0- 4810-0 0 OF10-0 8010-0 ZTEIO-0 9SI0-0 7610-0 €770-0 L 
7Z00-0 8£00-0 1200-0 9600-0 LOIO-0 8FI0-0 FL00-0 8010-0 F610-0 9170-0 6F£0- 0| 2070-0 0 9S00-0 8400-0 FOIO-0 9FI0-0 FIZ0-0 9 
0810-0 ZS10-0 PEIO-O SIIO-0 8800-0 0900-0 8Lh0-0 86£0-0 7LZ0-0 FOZO-0 LEIO-O TI80- 9090-0 8770-0 0670-0 9610-0 
8E10-0 EF1IO-0 IPIO-0 SII0-0 LL00-0 8£00-0 6L£0-0 Z8E0-0 ZSZ70-0 9LIO-0 8890-0 £120-0|09S0-0 P7LO-0 8E70-0 v 
6010-0 9€10-0 FS10-0 8E10-0 ZO10-0 1900-0 6670-0 FZE0-0 09E0-0 86Z0-0 8ZZO-0 0790-0 TZLEO-9 7670-0 € 
¥800-0 9600-0 O£10-0 #S10-0 ZS10-0 9Z10-0 6800-0 OFZ0-0 F870-0 ZEEO-0 9870-0 O8S0-0 0190-0! 6160-0 
0900 - 0 9400: 0 0 vrl0-0 8F10-0 vrlo- 0 0 €810-( 0 0 ) 8820- 0 0 Ist0-0 creo: 0 8ct0- 0 0 0890: css0- 0 ‘0 I 

| (syapou 

8I Ll 91 SI tl Ol 6 8 L 9 14 
sjuiod 
uo sjulog | 
| 


L ANY 9 sTaGOW 


ATaVL 


AONANTANI NOLLOIWAG LSAL 


9 1300W 


| 
| 
| 
| | 
| 
| 


WING DEFORMATION 


LOW ASPECT RATIO 


611-0 660-0 160-0 
OFI-O IST-O L8I-0 £07-0 
090-0 
607-0 


0070: 
9900 - 
L6t0- 
Pico: 
8870: 
81Z0- 
7790: 
6950: 


180 
LLI 


DS 


790 
£91 
897 


‘0 190- 


LLP: 
‘0 $SO- 


lel 


‘0 
8SZ- 


0 687 


0 


8tt-0 99¢ 


0 $90-0 290-0 OFI 


‘0 OO1-0 


1L0-0 L6P 
LS1-0 888- 


0 £9¢-0 


L971 -0 
7660-0 
8bL0-0 
9690-0 
0791-0 
+160-0 


ol 


Ol 


julod 108 


‘0 LOZT-0 E81-0 LHl-O 66E-0 £9E-0 OFE-O $67-0 


68t-0 
I 961- 
‘0 LOF-0 
0 8tZ-0 
691-0 


6 


NOLLIGNOOD 


SOOT 
0890-0 
81S0-0 
PSSl-O 
Strl-0 
OLT1-0 
OLTL-O 
+160-0 


6 


6710-0 
7L90-0 
OLIT-O 
ttOl-O 
P880-0 


8 


OSE-O $67 
£00-1 
8h1-0 OSI 
8tt-O0 LLZ7-0 607 
709-0 966-0 S9E 
$9L-0 092-0 OLL 


8 L 9 


IMIG 


A 


ZILO-0 
9011-0 
Z901-0 
9P01-0 
6001-0 
760-0 


0890-0 
0180-0 
PO80-0 
6£80-0 
1880-0 
1680-0 


L 9 


-0 6£5-0 £95: 
‘0 L6L: 
-0 


16: 


687: 
I 


v 


0 Sss-O Lts 
1 9LL-1 679- 
$7Z-0 
0 $69-0 709-0 IIS 
801-1 


TS81-0 


98S1-0 SSsSl-0 66h1-0 


(Speo] OOT O} BNP ‘saysul ul UMOYS SON|PA) 


Al 


IONANTANI NOLLOI144G0 


‘0 
66¢°1 
0 


KAM 


juiod 


i 

| 

’ 

_ 4 

Toe 

= 

| 
| | 

= j 

= 

a 
| 
NNN 
coocc 
| 
tT = 
| 
} 
= 
a) 
f 
3 
H 
| 
2 
| 
| i 
| 
| 
| 
| 


NODE 


0:28 
NODE 3 [2 
a ro) 
og \ € 
4 S I 
6 
Za 24 
2 
oa 
ad 
0-08} & 
ra) 
2 14/ 715 
0-04 | 
0 SUPPORTED 0 
(2) MODEL 4 (b) MODEL 4 
f Ficure 4. Deflected surface shape. Uniform load. Models 1 and 4. 
22r 
20F0 
é 
ad 
FiGurE 5. Deflected surface shape. Uniform load. Models 1, 2 and 5. 
NODE 
5 
MODEL 2 
ra 
Sw oY NODE 
1-0} 82 12 MODEL 
4 0-44 MODEL { 
$3 
O6+ Ea O3- Od \ \ \ 
ra) \ \ \ 
4 + 4 \ 
Ficure 6. Deflected surface shape. FIGURE 7. Deflected surface shape. 


Uniform load. Models 2 and 3. Uniform load. Models 6 and 7. 


LOW ASPECT RATIO WING DEFORMATION 


here are understood to be directed normal to the wing mid-plane.) To allow for an 
efficient presentation of data, and since the analytical methods to be evaluated 
necessarily produce a symmetric set of displacement influence coefficients, the 
experimentally derived cross-coefficients have been averaged to yield a symmetric 
set. The results are shown in Tables I-IV as displacements in inches, due to 100 Ib. 
loads. For economy in presentation the results for Models 1 and 2, 3 and 4, and 6 
and 7 have been paired. In each case, advantage has been taken of symmetry by 
tabulating only the values on and to one side of the main diagonal. 


Further operations on the test data are useful in the interests of providing 
concise parameters for the test-theory comparisons. There are numerous such 
parameters that can be obtained from the complete set of displacement influence 
coefficients and which represent certain basic properties of the elastic system. One 
is the fundamental wing vibrational frequency; another is the deformed shape under 
distributed loads. Since this study had an ultimate concern with static 
aeroelasticity, the latter parameter was selected. 


In particular, the chosen distributed loading case involved the simultaneous 
application of 100 Ib. loads at all spar-rib intersection points. The displacement at 
each reference point on each wing for this loading case, henceforth referred to as 
the uniform load condition, is listed in Table V. Graphical representations are also 
shown (Figs. 4-7) in the form of carpet plots of the deflected surface shape. Each 
plot typifies a geometric characteristic. 


Figure 4 illustrates the results for Models | and 4, the rectangular plan form 
wings. Results for the family of swept wings. Models 1, 2 and 5 (0°, 45°, and 30° 
sweep, respectively) are shown on Fig. 5. A comparison between Models 2 and 3, 
where the skin thickness is doubled, is given in Fig. 6, and Fig. 7 compares the delta 
plan forms. Models 6 and 7. It is noteworthy that the deformed shapes of Models 6 
and 7 differ significantly in consequence of the difference between internal member 
orientation. Not all of the node point locations on Models 6 and 7 exactly 
correspond, but this factor should not heavily influence the level of agreement in 
the uniform load comparisons. 
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The Static Aeroelastic Deformation of 
Slender Configurations* 


PART II: Some Calculations on Slender Plate Aircraft, 
Including Non-Linear Aerodynamics 


G. J. HANCOCK 


(Department of Aeronautical Engineering, Queen Mary College, 
University of London) 


SUMMARY: Calculations have been applied to a series of flexible plates (delta, 

gothic and ogee) to estimate their static aeroelastic behaviour. Assuming linear 

aerodynamics, the general tendencies of Part I are particularised. Assuming 

non-linear aerodynamics, it is shown that in general two positions of 

equilibrium are apparently possible at each speed. The stability of the models 
is not considered. 


1. Introduction 


In Part I a qualitative idea of the primary type of aeroelastic deformation of a 
slender configuration under linear aerodynamic loading, involving a lengthwise 
deformation only, is given for the trimmed cruising condition in level flight at 
various forward speeds. It is shown that the aeroelastic deformation may become 
very large at a finite speed when the elevator control force to trim the aircraft also 

ecomes very large. In this note some detailed calculations are carried out to 
supplement these statements with quantitative results, and to extend these ideas to 
cover cunfigurations under non-linear aerodynamic loadings (i.e. including leading- 
edge vortices). 


The procedure followed here is now indicated. 


(i) The equation for the lengthwise deformation is derived for general 
distributions of weight and stiffness, including expressions to represent the non- 
linear aerodynamic loadings associated with slender configurations, based on the 
assumption that the whole airframe can be treated as a free-free beam. This 
assumption may not be a severe one for a slender configuration, which will comprise 
a centre body fuselage supported by wings of a slender plan form. The main 
difficulty is introduced by the present lack of knowledge regarding the non-linear 
aerodynamic loading on streamwise cambered slender wings. A linear theory” is 


*A summary of part of this work was read at the meeting of the AGARD Stability and 
Control Panel in Brussels, April 1961. 


Received January 1961. 


372 The Aeronautical Quarterly 


a 
e 
t 
i 
\ 


av: 

vo 

wi 

the 

the 

be 

| be 

He 

na 

th 

in 

th 

ae 

ec 

Sc 

fc 

sé 

re) 

SI 

d 

a 


STATIC AEROELASTICITY 


available for this case, but this theory is only applicable when the leading-edge 
vortices are absent. Non-linear theories have been established for delta wings, 
without camber, to illustrate the effects of these leading-edge vortices. In this paper, 
the non-linear form of the loading is retained in general form (in the derivation of 
the fundamental deformation equation), so that, as more aerodynamic information 
becomes available, substitution of more exact theories may be easily accomplished. 
being of an arithmetical nature, without altering the approach to the problem. 
However, for the purposes of calculation here an intuitive representation for the 
non-linear aerodynamic loading for general shapes with camber is made, based on 
the existing information referred to so that, although the quantitative results derived 
in this paper are dependent on this assumption, it is hoped that these results reflect 
the general aeroelastic behaviour of slender configurations. 


(ti) The general equation derived in (i) is solved exactly for the case of a 
trimmed delta plate of constant thickness with both the linear and non-linear 
aerodynamic loading, for both its deformation and control force characteristics. 


(iii) A simple numerical technique for an approximate solution to the general 
equation derived in (i) is suggested, and is applied to the trimmed delta plate 
solution derived in (ii) to check the accuracy of the method. 


(iv) The numerical technique, introduced in (iii), is applied to the cases of 
gothic and ogee trimmed plates with linear aerodynamic loading. 


(v) The gothic trimmed plate with non-linear aerodynamic loading (i.e. in the 
presence of leading-edge vortices) is investigated. 


It is seen that the actual calculations are restricted to plates of different plan 
forms with constant thickness. Although it may be argued that this is not repre- 
sentative of a typical slender configuration in regard to either weight distribution 
or stiffness distribution, at least the important parameters are inherent to these 
simple cases, which are consistent in their relation between weight and stiffness 
distributions, particularly in the respect of the effects of change of plan form. In 
any case the numerical method, mentioned in (iii), is applicable to the more 
general configuration. 


Finally, at this stage it is worth noting some of the implications of the 
applications of linear and non-linear aerodynamic theories. In the interests of 
efficiency the slender aircraft may cruise without any appreciable leading-edge 
vorticity, and this can be achieved by a design incorporating a spanwise camber“. 
So the application of linear aerodynamic theory at cruising speeds for such shapes 
will be valid. It will still be essential to investigate the deformation characteristics 
below the cruise condition when the leading-edge vortices will play an increasingly 
important role. However, this aspect is not covered directly here, since we shall 
consider a given plate without any spanwise camber so that the leading-edge 
vorticity will always be present when lift is generated. 
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y 
© | 
| — —— 2s (oc) 
Lx) LOAD 
(ac 
(2) WEIGHT 
DISTRIBUTION 
Ficure 1. Slender aircraft notation. 
NOTATION 
° wing plan form 
& Pp 


I(x) load distribution 
1x), L(x) generalised functions in the load distribution (equation (8)) 
K=tany 
P trailing edge control force 
w(x) weight distribution 
W all-up weight 
V relative free stream velocity 


V. maximum trim speed 


z=( (x), wing camber shape 
E=x/c 
a(x)= —{ (x), incidence distribution 
A speed parameter (equation (30)) 
A, weight parameter (equation (31)) 
o weight distribution (equation (31)) 
y semi-apex angle of the plate 


Poisson’s ratio 
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2. The General Deformation Equation 


The slender configuration is shown in Fig. 1. The origin of a right-handed 
Cartesian set of co-ordinates is taken as the apex of the wing in the zero-lift 
position. The x-axis is in the relative stream direction and the y-axis denotes the 
spanwise direction. The relative free stream velocity is denoted as V. The plan 
form of the wing on the plane z=0 is given by 


It is assumed that s’ (x) > 0 for all the plan forms considered. The attitude of the 
wing to the free stream is denoted by 


which is assumed to be a function of x only, since in this paper interest is restricted 
to the longitudinal bending of the aircraft. The weight distribution is denoted by 
o (x, y) Ib./ft.?, so that 


+s OO 


weight 
w= init length | (x, y) dy (3) 
(x) 
and so W= total weight w(x)dx. . (4) 


The effective longitudinal bending stiffness distribution may be represented as E/ (x), 
which can be calculated for any given structure. 


The aerodynamic loading is now considered. The primary interest is the lift 
distribution along the chord. According to linear slender-wing theory”, assuming 
that the configuration is sufficiently thin for the theory to be applied, the lift 
distribution / (x) Ib./ft. along the chord is given by 


dx 

It is seen that, on integration of equation (5), the total lift is given by the conditions 

at the trailing edge only. In the presence of leading-edge vortices, equation (5) is 

seriously inaccurate, as demonstrated on a delta flat plate (i.e. symmetrical sharp 

edged section) wing’. The formula 


L(x)= pV? + 48 0) [VO] (6) 


is reasonably accurate for this type of wing, i.e. when s’ (x) and Y (x) are constant, 
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On the basis of both theoretical and experimental evidence. No definite information 
is available at present to cover cases where s’ (x) is not constant (i.e. curved leading 
edges) or where (’ (x) is not constant (i.e. chordwise cambered wings) so that some 
assumptions must be made in order to proceed with the general distortion problem. 


Now equation (5) is completely general, with regard to s’ (x) and ¥ (x), but is 
limited by its linearity. Equation (6) is a reasonable non-linear equation when s’ (x) 
and ( (x) are constant, and the first term (or the linear term) of equation (6) is then 
identical to equation (5). Therefore it is suggested that an approximation for the 
general non-linear case could be written 


L(x) pV? (x) (0x) + (2) ds| (7) 
a. 


which includes equation (5) for the linear first term and satisfies equation (6) when 
s (x) and ( (x) are constant. It is admitted that there are several formulae which 
could be used to represent the non-linear term in equation (7), applying similar 
arguments to those given already. However, equation (7) satisfies one most 
necessary criterion for general usage, namely, that the load remains finite at the 
trailing edge if s’ (x) tends to zero there. 

Since the accuracy of equation (7) is an unknown quantity, it will be assumed 
that the form of the aerodynamic loading may be written 


d 


l(x)= —pV* 


Eo U(x) + 1.00) YO) | : (8) 


and the fundamental deformation equation will then be derived in terms of /,(x) and 
L(x). As more information becomes available, accurate expressions may be 
substituted for /,(x) and /,(x) but, for the purposes of calculation in this paper, it 
will be assumed that 


x 


and 1, (x)=4 | s(x) dx, 
by comparison of equations (7) and (8). 


It is assumed that the aircraft is trimmed for level flight by a trailing edge 
force P, which is induced by the inclination of the elevators. 


In level flight the overall conditions that 
lift = weight ‘ (11) 
must hold. 
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The first of these equations states that 


from equations (8) and (4), remembering that /,(0)=/,(0)=0. The second 
condition is 


: 
py iy L(x) (0) + LQ) w(x) {c—x}dx=0, . (13) 


by taking moments about the trailing edge, where w (x) is defined in equation (3). 


As explained in the Introduction, it is assumed that the structure may be treated 
as a beam, so the deformation equation is 


| EI’ (x)) + pV? 


J 


together with the boundary conditions 


EIt”|,..=0 
EI” |,<. =0 (16) 
dx 
ax 


Equations (15) and (16) state that no bending moments are applied at the apex or 
the trailing edge. Equation (17) states that no shear force is exerted at the apex. 
Equation (18) states that the control force P acts at the trailing edge. 


On integration, equation (14) gives 
+ pV? + LQ) | w Gd¥=0. . (19) 
a. 
Equation (19) satisfies the boundary equation (17), since the loading is zero at the 
leading apex. Equation (18) is also satisfied if equation (12) holds. 


On further integration, equation (19) gives 


+ pV? 11,00 + dxt | dx w(x)dx-0. . (20) 
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Equation (20) satisfies the boundary condition (15). Using equation (20) it can be 
shown that boundary condition (16) is equivalent to equation (13). 


If we now introduce the non-dimensional chordwise co-ordinate 


E=x/c (ie. OS EX 1) 

then equation (20) becomes 
Elz‘(é) pVic? | LEO) 2(O dé+e* | dé, w dé, (23) 


0 0 0 


where the prime (") now denotes differentiation with respect to €. The remaining 
boundary condition (i.e. equation (16)) to be satisfied by the solution of 
equation (23) is 


assuming that E/(1)> 0. Thus an integral-differential equation has been obtained 


for the incidence distribution 2(£) of the flexible shape in terms of the known 
distributions E/ (£), w (£), 1,(€) and 1,(€). 


3. Slender Delta Plate Wing (Linear Aerodynamics) 


We now consider the simple case of a slender delta plate wing, of constant 
thickness t and uniform weight distribution «, trimmed for level flight by a control 
force P at the trailing edge, as explained in Section 2. In this case, from equations 
(1), (3), (9), and (10), 


W (x)= 208 (x) =20Kc€ : : (26) 


2Et’s (x) 


EI (x)= 


TS°(x) = | 


L(x)=4 | s(x) dx= | 


where K is a constant and equal to tany, where y is the semi-apex angle of 
the plate. 
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Substitution of equations (25), (26), (27) and (28) into equation (23) gives 


Al (2+ “3 (29) 
where A and A, are constants given by 
and 
The boundary condition for the differential equation (29) is, by equation (24), 


First the linear problem only is considered, so equation (29) becomes, on 
ignoring the non-linear terms, 


| Sade+A,2. (33) 

It is seen immediately that one solution is 
~constant= (34 


This is due physically to the fact that the distribution of aerodynamic loading 
balances the weight distribution, so that no deformation exists. Also the trailing edge 
control force P is zero (from equation (i8)). This problem is therefore similar to 


(4) APPROXIMAT 
( EXACT , 
a (SECTION 3) 
Z 
0-2 04 0-6 0-8 1-0 


FIGURE 2. Deformation mode of delta plate (linear aerodynamics) at the maximum trim speed. 
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Figure 3. Two modes of equilibrium of delta plate (linear aerodynamics) 
at the maximum trim speed. 


Model 4 discussed in Part I of this paper”, in which it is also indicated that a speed 
may exist at which there are an infinite number of positions of equilibrium. Thus 
the solution given in equation (34) is not necessarily unique for all speeds. 


The general solution of equation (33) may be obtained by the method of series. 
and it may be shown that 


(£) z,—(Aa,—A,) (35 (35) 


where z, is a constant and represents the incidence of the nose. The boundary 
condition, equation (32), gives 


1 <A A? 
(Aa,—A,) (3 - + ) =0. (36) 


Therefore, either a,=Ay/A, 


which is the solution given by equation (34), or 
=0. . (37) 


When condition (37) is satisfied, z, is indeterminate, but the control force P is non- 
zero and depends on z,. The value of A which satisfies equation (37) is 


which gives, by reference to equation (30), a critical speed at which the deformation 
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amplitude is indeterminate. The angular mode of the deformation 2, (¢), relative to 
the undistorted shape at this critical speed, is given by 


g 336 2° 


and this function is shown in Fig. 2. The amplitude of this mode is governed by 
the magnitude of the trailing edge force P by reference to equations (12) and (35). 
If P< W (the all-up weight), then z(1)>0. If P> W, then z(1) <0, but this 
condition on the control force is impractical since it implies that all the lift is 
generated on the controls. However, for a sample calculation, for the sake of 
simplicity, we may take the case when P=W; then (1) is zero, so that, from 
equation (35), 

KrpV,? 

where V, is the critical speed. This equilibrium position is shown in Fig. 3, together 
with the attitude when P is zero. 


A more concise solution of equation (33) may be obtained directly in terms of 
Bessel functions; this is shown in Appendix I. 


4. Slender Delta Plate Wing (Non-Linear Aerodynamics) 


The case with the non-linear aerodynamic loading is now considered, that is, 
the full equation (29). Experience in the solution of the linear case (equation (33)) 
shows immediately that the non-linear case also has a similar trivial solution, that is, 


2 
x | ] + x constant = 


(39) 
nK A KzpV° 


There is only one positive attitude z, which satisfies equation (39). Physically there 
is again no deformation because the aerodynamic loading and weight distribution 
balance. The presence of the vortices does not alter the distribution of loading; 
only the intensity of loading is affected, so that the plate with the leading-edge 
vortices takes up a smaller incidence than the plate without leading-edge vortices. 


The determination of an equivalent to the critical speed, including the effects of 
the leading-edge vortices, by the general solution of equation (29) is difficult and 
has not been fully solved. The difficulty arises because of the (z| term, when 2 (£) 
might be expected to change from a positive value at €=0 to a negative value at 
€=1, for large values of the control force P. However, if it is assumed that z does 
not change sign, the method of series can be applied and the solution will then only 
be valid when the control force P is less than W, which is not a severe practical 
imposition. 
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The series solution of equation (29), assuming that 2 (£) is positive, is 


2 A 4 


16 


+ OK °° * 


3 


and the boundary condition gives 


A? 
Obviously equation (39) represents one solution of equation (41), but equation (41) 
opens up many further interesting points. For example, it is now possible that for 
each value of A there are, at least, two values of 2, which satisfy equation (41), 
obtained by equating both brackets separately to zero. 


Now there is at least one value of A when both brackets are identically zero; 
this is given by the simultaneous conditions that 


(1 2) A 


n 9 


By reference to equation (37), and its solution in equation (38), these two conditions 
(42) become 
(1 


(43) 
and A + 33 | 


Thus, by solving these two equations, we obtain both the incidence 2 (=z,) and the 
speed condition .\ at which the trailing edge boundary condition is “doubly” 
satisfied. 


As an example, if we take A, =3-6 and 2/(*K)=2°5, then the unique solution 
of equations (43) is 


(44) 
a(é)=2,=0:137 | 
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When A < 19-6, two equilibrium positions may exist. If we take the speed given 
by A= 16'S, the first bracket in equation (41) gives 


2, (1+2-5z,)=0-22 


with its associated deformation solution 


The second bracket in equation (41) gives 


\! 36 + g +202, + 502, 0. 
The solution of this equation is approximately 
. . (46) 


so that, on substitution into equation (40), a second deformation solution is 
z (€)=0-20 [1 — 0-75 + 0:69¢* +... .. (47) 


This series is convergent and z(£) > 0 for 0 € <1, so that the solution is valid. 
Thus, equations (45) and (47) are the two equilibrium positions for the speed given 
by A=16°S. 

It is also apparent that in a case where A > 19-1 two finite solutions again may 
exist. For example if A =22 then the two solutions are 


a (£)=0-°125 | 
and, approximately, 2(€)=0°10 [1 + 0-97& . 


Again, the series for the second solution is convergent and positive for 0<= ¢ = 
so this solution is valid. 

Thus it has been shown that a delta plate has two positions of equilibrium, at 
least over part of its speed range, when non-linear aerodynamics are assumed. The 
extreme values of A, either very small or very large, have not been considered, since 
their solution is not amenable to the technique applied here. It is established, 
however, that no reasonable speed exists at which an infinite number of equilibrium 
positions are possible, depending on the magnitude of the control force, as with 
linear aerodynamics. 

The question which is now posed is whether both positions of equilibrium at 
any speed are stable. This can only be answered by reference to stability, which is 
outside the scope of the present paper. 
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5. A Numerical Method 
We consider first the linear form of equation (23) which is 


| 2() dé +e | de | (49) 

Now, in principle, if the known functions E/ (£), 1,f€) and w (¢) can be expanded as 
polynomials in €, then a solution can be obtained by the method of a series 
expansion. Some initial difficulty occurs if the polynomial for [E/(¢)]~' has a 
radius of convergence less than unity, for then a transformation of variable is 
necessary. This occurs, for example, on an ogee plate (see Section 7). However, the 
convergence of the series for z(€) is extremely slow when £ = | (where the boundary 
condition has to be applied), so the method becomes tedious. In a practical case it 
is possible that only discrete values of E/ (€) and w (€) at particular points along the 
chord are available, so that the application of the series method becomes even more 
formidable. It is therefore desirable that a simple numerical technique should be 
introduced. Fortunately equation (49) is expressed in a suitable way for a simple 
step-by-step process. Starting at €=0 with a value of 2 (0) then 


€ € 


EI (2) | | dé, w(é)dé, . SB) 


where < is a small value. Thus 2’ (<) can be determined, since the other terms are 
calculable, and then it follows that 


The next step is the calculation of 2’ (2) from the equation 


dé,|w(Odé.. (52) 


€ 


This process can then be carried on until z’ (1) is reached. This numerical method 
is easy to apply. The functions 


(n+1)€ neé 
EI (n:), and | dé, | dé 


ne 0 0 


are predetermined from the given distributions for all the points nz, where 
0 =n<=N and Ne=1. In the author’s applications it is assumed that ¢=0-1 and, 
as will be seen, this seems to give reasonable accuracy. The method described in 
equations (50), (51) and (52) can be applied rapidly, taking about 15 minutes to 
reach €=1 on a desk machine. 
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The basis of the method is the determination of the initial value 2(0) which 
ultimately satisfies boundary condition 2z’(1)=0 after the application of the 
numerical technique. One could repeat the process for various values of 2 (0) and 
so find the correct z (0) which gives z’ (1)=0. However, since equation (49) is linear, 
the problem may be reduced to two independent problems, by noting that the 
complementary solution of the equation 


6 
may be added to any particular solution of equation (49). Thus, the numerical 
technique can be applied to equation (53) to find a function z. (€), assuming initially 
a, (0)=1 . (54) 
and so determining z.’(1). Then the technique can be applied to equation (49) to 
find a function assuming 


. . . (55) 


and so determining z,'(1). Therefore, if the constant a, is determined from the 
condition 


. . . (56) 


then the solution of equation (49) which satisfies the boundary condition 2‘(1)=0 is 


One important point which arises directly from equation (56) is that, when 


a, becomes infinite, assuming that z,'(1) is non-zero. This, therefore, is the con- 
dition for the maximum trim speed at which a total loss of control power occurs. 
The deformation will tend to infinity as the speed approaches this critical condition. 
This is the phenomenon discussed in relation to Model 5 of Part I (Ref. 1). Thus 
the determination of the maximum trim speed in this linear case may be considered 
(equation (53)) independently of the overall deformation problem. Numerically the 
maximum trim speed is determined by applying the numerical method to equation 
(53) at increasing speeds until the boundary condition (equation (58)) is satisfied. 


Although the delta plate wing is a special case, it illustrates the points already 
discussed. First, the accuracy of the numerical method can be assessed by 
comparison with the exact solution given in Section 2. The numerical method gives 


A=32‘5 . ‘ (59) 


for the maximum trim speed. This is approximately 2 per cent different from the 
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exact value (equation (38)). A comparison of the approximate and exact deforma- 
tion modes is given in Fig. 2, indicating reasonable agreement. At speeds below the 
maximum trim speed, the application of the numerical method, as described in the 
discussion of equations (54), (55) and (56), gives the result 


Thus there is no deformation, which agrees with the exact solution. 


The numerical technique described in equations (50), (51) and (52) may be 
extended to the case with non-linear aerodynamics (equation (23)), although the 
simplifications described in equations (53)}(56) for the linear case cannot now be 
applied. Thus, at a particular speed an initial value of z(0) has to be guessed and 
@’(1) calculated. By the variation of the values 2 (0), the correct value of (0), 
together with its associated mode 2 (£), can be found which satisfies the boundary 
condition 2‘(1)=0. A sample calculation is indicated in Appendix II. 


6. Gothic Wing (Linear Aerodynamics) 
The next case which is considered is that of a plate wing with a gothic plan 
form (Fig. 4) given by 


s(x)=cKé (1 5) . 
In this case the differential equation is 
¢ 4 ¢2 é ¢ 3 
5) Afi (1-5) 2@de+a,4 (1 (62) 
where A and A, are defined in equations (30) and (31). The application of the 


numerical technique gives the maximum trim speed as 


A=44'5._. . (63) 


The equilibrium shapes of the wing at various values of A below the maximum trim 
speed are shown in Fig. 5. It is seen that, as A tends to the critical value, the 
deformation increases. The value of the control force P as a function of A is shown 
in Fig. 6 and its divergence is indicated. 


7. Ogee Wing (Linear Aerodynamics) 
To investigate further the effects of plan form on deformation characteristics, a 
plate wing of ogee plan form, trimmed for level flight, is now considered. 


The plan form (Fig 7) is given by 
. ‘ . (64) 


This introduces a slight modification to the plan form of the delta plate investigatec 
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in Section 3. For this case the differential equation is 


0 


with the usual definitions of A and A,. As mentioned previously, this equation 
cannot be solved as it stands by the method of series, since the term (1 + ¢—£*)~? 
has a radius of convergence less than unity. However, the numerical technique can 
be applied; this gives the maximum trim speed as 


The equilibrium shapes of the wing at various values of A below the critical values 
are shown in Fig. 8. The variation of the control force P as a function of A is 
shown in Fig. 9. It is seen that the rate of divergence of the control force P is rapid 
in the region of the critical condition. The loss of control effectiveness for trim has 
more significance for this ogee wing than for the gothic wing, in Section 6, where 
control would be lost in practice at speeds far below the critical. 


8. Gothic Wing (Non-Linear Aerodynamics) 


The case of the gothic wing considered in Section 5 is now extended to 
investigate the effects of assuming non-linear aerodynamics. The differential 
equation is, by reference to equations (62), (23) and (10), 


| (67) 


It is seen from equation (67) that the shape of the wing must be further defined 
before any calculations may be carried out. It is now assumed that K=4. The 
procedure, as explained previously, is the calculation of the equilibrium position of 
the wing for low values of A, repeating the calculations for increasing values of A. 
This has been applied to equation (67) and the two important points which arise 
from these calculations are as follows :— 


(i) For low values of A, two positions of equilibrium are found. 
This is similar to the delta wing (Section 4). 


(ii) Ata certain value of A there is only one position of equilibrium 
but, above this speed, no position of equilibrium is possible. 
This is essentially different from the de'ta plate case, where 
positions of equilibrium could always be found. There is really 
no resemblance to the case of the gothic wing with linear aero- 
dynamic loading, where the deformation becomes infinite at a 
critical speed and there is divergence of the control force. 
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However, in this case with the non-linear loading, there is no 
divergence of the control force and the deformation is always 
finite but, above a certain speed, the aircraft cannot be trimmed. 


The limiting speed for equation (67) is 


The limiting distortion is shown in Fig. 10, together with the two positions of 
equilibrium at A = 18. 


The limiting speed with linear aerodynamics is A = 44:5 (equation (63)), so the 
limiting velocity is reduced by approximately 50 per cent by the introduction of 
non-linear aerodynamics. This may be significant for the aircraft which is designed 
to cruise without any appreciable leading-edge vorticity but will have to fly below 
its cruise conditions when the effects of leading-edge separations become more 
important. 


9. Discussion 

Deformation characteristics of slender plates have been discussed, assuming 
both linear and non-linear aerodynamics. Only the implications of these loadings 
have been discussed and there is no emphasis on which loading is the more practical. 
This is a matter for further work, but at least it is hoped that the difference between 
these loadings has been satisfactorily established. 


With the assumption of linear aerodynamics a maximum trim condition occurs, 
with the aeroelastic deformation and the control force tending to infinity at the 
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Ficure 10. Equilibrium shapes for gothic plate (non-linear aerodynamics). 
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critical speed. The rate of loss of control effectiveness to trim at speeds approaching 
the maximum trim speed depends on the out-of-balance between the static inertia 
and aerodynamic loads. In the limit, when these two load distributions exactly 
balance, as for the delta plate, there is no deformation below a critical speed, but an 
infinite number of positions of equilibrium are possible at this critical speed, 
depending on the magnitude of the control force. 


With the assumption of non-linear aerodynamics, a total loss of control 
effectiveness to trim apparently does not occur. Only two plate aerofoils have been 
investigated, a general delta (i.e. with variable aspect ratio) and a particular gothic 
(aspect ratio 3/4). Both plates have a “critical” speed at which there is only one 
equilibrium attitude, which is finite, and both plates have two positions of 
equilibrium at any speed below this “critical” speed. However, the delta plate also 
has two positions of equilibrium for each speed above the “critical” speed, while 
the gothic plate cannot be trimmed (i.e. no equilibrium state is possible) above its 
“critical” speed. Further work is needed to establish the relationship between the 
weight distribution and the aerodynamic loading when an upper limit on speed 
exists due to deformation. The study of the simple models used in this paper has 
enabled us to determine the trends and the fundamental aspects of the effects of the 
static aeroelasticity of slender configurations in level flight. The incentive to extend 
these ideas must come from the need to know the implications of more practical 
distributions of weight, bending stiffness, and so on, than are assumed in this note. 
The numerical method (Section 5) covers these more complicated distributions and 
their solution is a matter of routine calculation. 


The significance of all these results is bound up with the relationship of static 
aeroelastic and dynamic aeroelastic effects as defined in Part I. We have only been 
concerned with the steady manoeuvre of trimmed flight of a slender configuration in 
this note. To complete our understanding of the overall picture, the distinct problem 
of the stability of this equilibrium solution is at present under consideration. 


Appendix I 


The solution of the equation (33), that is, 


may be expressed directly in terms of Bessel functions. 
On differentiation with respect to ¢, 
fa” ‘ (70) 
where the boundary conditions are 
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The substitution E=u?/3 ‘ ‘ (72) 


in equation (70) gives 


where a,=4A/9, a,=4A,/9, and the prime (’) denotes differentiation with respect to wu. 
The solution of (73) which is finite for 0 <u <1 is 
a=A,J, {ur (a,)} +a,/a, (74) 
where A, is a constant. 


The boundary conditions are 


=0 where u=0 and u=1. (75) 
du 
dz 
Now =—/(a,)A,J, {u /(a,)}. 
w 
Therefore dz/du=0 as u—>0, since J,—>0. And, since dz/du=0 at u=1, then 
either A,=0 or J,{4/(a,)}=0. The solution A,=0 in equation (74) gives the equivalent A 
shape with no deformation. The fundamental solution of J,{/(a,)}=0 is /(a,)=3-83, - 
which leads to A=33. This agrees with the series solution (equation (38)). It is 
noted that higher mode solutions are available at approximately /(a,)=7-0, 10-2, ” 
13-3, ote. 
Appendix II 
APPLICATION OF THE NUMERICAL TECHNIQUE TO THE GOTHIC PLATE WITH ’ 
Non-LINEAR AERODYNAMICS 
Zz 
The equation to be solved is equation (67), namely 
with the boundary condition 3 


2’ (1)=0. 


The numerical technique which is to be applied to this equation is described in 
Section 5. This is a step-by-step process, from £=0 to ¢=1 and in steps of :=0-1. 
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TABLE I 
ne 
é €(1—£/2 | &(-€/4) €2 (1—£/2)dé 
(n- De 

0-0 | 0-000 1-273236 0-0 
01 0-095 1-363762 0000975 0-0185 
0-2 0:180 1-467103 0:0076 0°1177 
0:3 0:255 1-586038 0024975 0:2845 
0-4 0:320 1-724175 0-0576 0:4876 
0-5 1:886274 0109375 0-7011 
06 0-420 2:078753 0: 1836 0-9066 
07 0-455 | 2310409 0282975 10954 
0-8 0-480 2:593628 04096 12710 
0-9 0-495 2946330 | 0:564975 1-4504 
1-0 0-500 3-395296 0:75 16672 


To start with it is necessary to compute the functions 


ne 
4 (1—£/3) ag & - 


n~—De 


where 0 <= n< 10. These are indicated in Table I. 


A sample calculation for A==18 and .\,=3-6, assuming the initial value z(0)=0-4, is 
as follows :— 


1 18 
=()-4: 4 ): + 1-2 (0- 5)} = —0-022955 
2(0) =0-4; 2’ (0-1) (0011169) + 1-2 ( 000975) 0-022955 
0-1)=0-397705; 2’ (0-2)= 
0-18 [—0-3 (0-083367)+1-2 (0-0076)] =—0-088273 


a (0:2)=0-378878; 2’ (0-3)= [—0-3 (0:251073)+ 1-2 (0-024945)]= —0-177851 


z (0-3)=0-361093; 2’ (0-4)= [—0-3 (0-527978) + 1-2 (0-:0576)] = —0-278979 


32 
z (0-4)=0-333195; a’ (0-5)= [ —0-3 (0-895782) + 1-2 (0-109375)] = —0- 366625 


z (0:5)=0-296533; ¥ 0-)= {—0-3 (1-314990) + 1-2 (0-1836)] =—0-414707 


(0:6)=0-255052; 2’ (0-7)= 0 {[—0-3 (1:742001) + 1-2 (0-282975)]= —0-402594 


z (0-7)=0-214793; a (0°8)= [—90-3 (2-150483) + 1-2 (0-4096)] =—0-320051 
l 

z (0-8)=0-182783; (0°9)= {[—0-3 (2:541274)+ 1-2 (0-564975)]= —0-170529 


z (0-9)=0- 165730; [—0-3 (2:952497) + 1-2 (0-75)] = +0-285. 
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This process can be written down as shown, with the aid of a desk machine, in 


about 15 minutes. By varying = (0), the associated value of a’ (1) can be found. By 
plotting a (0) against a’ (1), the values of 2 (0) for which a’ (1) is zero can be determined, 
and so can the corresponding mode. 


REFERENCES 


Hancock, G. J. Static Aeroelastic Distortion of Slender Configurations. Part I. Aero- 
nautical Quarterly, Vol. XII, p 293, August 1961. 


Jones, R. T. Properties of Low Aspect Ratio Pointed Wings at Speeds Below and 
Above the Speed of Sound. N.A.C.A. Report 835, 1946. 


MANGLER, L. W. and Situ, J. H. A Theory of Flow past a Slender Delta Wing with 
Leading Edge Separation. Proc. Roy. Soc. A, 251, May 1959. 


MASKELL, E. C. and Weper, J. On the Aerodynamic Design of Slender Wings. Journal 
of the Royal Aeronautical Society, December 1959. 


394 The Aeronautical Quarterly 


Si 


| 
3. 
R 
4. 


The Correlation of Unstable Crack Length 
Data for Sheet Materials 


E. R. WELBOURNE 


(Technical Department, Royal Aeronautical Society) 


SUMMARY: A semi-empirical analysis of unstable crack propagation in thin 
sheets is given. An effective value of the crack tip strain concentration factor 
is used as the criterion for instability. Good correlation is achieved with 
experimental results from unreinforced sheets of aluminium and titanium 
alloys and high-strength steels, The interdependence of geometry and stress- 
strain curve properties in determining the unstable crack propagation 
characteristics of sheet materials is illustrated in an appendix. 


1. Introduction 


Fracture mechanics has been a subject of some interest to the aircraft structures 
engineer for a number of years past. However, the problem of predicting the 
relationship between unstable crack length and stress for a fatigue crack in a 
structure is still far from solution in any general case. All theoretical analyses of 
the basic problem depend on at least one empirically determined parameter’: ? * 
and many purely empirical methods are in use. The mechanism of stable crack 
propagation is not properly understood although a tentative model has been 
proposed 


Deficiencies in the basic theory and reliance on purely empirical representation 
have tended to obscure certain important features of the experimental data which 
have been collected. This is particularly so in regard to unstable crack length results 
obtained from unreinforced sheets. In this note an attempt is made to show how 
such results can be correlated with a simple semi-empirical analysis, and to indicate 
the factors influencing the unstable crack propagation characteristics of various 
materials. 


2. Basic Theoretical Considerations 


In order to be quite clear about the initial assumptions of the analysis method 
proposed, brief consideration of basic principles is necessary. 


For unstable propagation of a crack in any material two conditions must be 
satisfied. 
(a) The strain at the crack tip must be high enough to cause fracture 
of the material. 
(b) The energy released during an incremental advance of the crack 
tip must just exceed the energy absorbed in this advance. 
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In a material which fails elastically, i.e. without experiencing plastic deforma- 
tion, these conditions are satisfied simultaneously; in a real material condition, 
(a) is satisfied before condition (5), in general, giving rise to a period of stable crack 
growth. However, if the stable crack growth period is neglected, then the two 
conditions are assumed to be satisfied simultaneously. By this assumption either 
one is then a sufficient condition for instability. Existing analyses invariably 
disregard the slow growth period and, in view of this, it is perhaps surprising that 
so much work has been based on theoretical analysis of condition (b). Approximate 
theoretical analysis of condition (a) is considerably simpler and less restricted in 
application. Furthermore, such an analysis requires only the normal engineering 
concepts of stress and strain, avoiding the obscurity of some of the material 
properties encountered in an energy analysis. For these reasons an approximate 
theoretical analysis of condition (a) is used in correlating the experimental results. 


NOTATION 


a_ length of sheet in direction of applied load (in.) 
b width of sheet (in.) 


t thickness of sheet (in.) 


modulus of elasticity of sheet material (Ib./in.’) 

f, tensile strength of sheet material (Ib. /in.*) 

f, 0-2 per cent proof strength of sheet material (Ib./in.*) 
Emax nominal elongation of sheet material 

material characteristic 

r effective radius of crack tip (in.) 

f average stress on net area of sheet (Ib. / in.) 

f, average stress on gross area of sheet (Ib. /in.*) 

€ average strain corresponding with stress f 

k effective strain concentration factor for crack 


X, unstable crack length (in.) 


3. Analysis 
3.1. Basic EQUATION 


Unstable propagation of the crack is assumed to occur when 


i.e. when the crack tip strain reaches the maximum elastic-plus-plastic strain of 
which the sheet material is capable. For convenience in calculation, maximum 
stress and maximum strain are taken to be coincident on the stress-strain curve. It 
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Ficure 1. 
is further assumed that nominal material properties can be applied when considering 
the highly localised region around the crack tip. This is possible since k values are 
determined empirically. 


3.2. DETERMINATION OF k 


In calculating values of the effective strain concentration factor k, two 
preliminary assumptions are made. 
(i) The effective radius of the fatigue crack r is a constant for 
the material. 


The expression for k found by elastic analysis is valid for large 
plastic strains. 


(ii) 


On this basis, for a central crack in a sheet of finite width subject to a uniform 
tensile stress normal to the line of the crack, 


(3; 
using the well-known expression for the infinite sheet modified for finite width, 
as in Ref. 5. Fig. 1 shows k as a function of x,/b for some arbitrary values of 
band r. For values of x,/5 less than about 0-4, equation (2) agrees closely with the 
strain concentration factor derived from actual tests on cracked sheets by using 
equation (1). For values of x,/5 greater than this, however, the decrease in k pre- 
dicted by elastic theory is not in agreement with k values derived from actual test 
results, which remain approximately constant, as indicated in the figure. As far as 
is known, no independent evidence of this effect exists and its interpretation in terms 
of a constant value of strain concentration is open to question. However, within 
the framework of this analysis, it can be accepted insofar as it must give the 


(x,/b)} (2) 


right answer. 
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FiGuRE 2. & for aluminium alloys. 


Figures 2, 3 and 4 give values of k for three materials computed on the basis 
given in the previous paragraph, using a value of r for each material chosen 
empirically to fit the available unstable crack length results. Some limitations on 
the application of these curves must be noted. 


No allowance is made for the variation of strength of centrally cracked plates 
with thickness, which arises from buckling of the crack edges. The curves given 
correspond with a sheet thickness of about 0-064 in. and predictions of unstable 
crack length in thinner sheets based on these k values may be slightly unconserva- 
tive. Conversely, conservative results will be obtained for thicker sheets. 


The effect of sheet aspect ratio has been neglected by assuming uniformity of 
the stress field to be achieved between the ends of the sheet and the crack region. 
This implies a value of a/b of about 2 or greater. Results obtained from square 
sheets tend to be influenced by the end conditions. 


Values of k are quoted for sheets up to 30 in. wide. Data for sheets wider than 
this are scarce for aluminium alloys and virtually non-existent for steel and 
titanium. Figures 3 and 4 are based on results from 10 in. wide sheets only. 
Extrapolation of the analysis to sheet widths over about 30 in. leads to substantially 
conservative results when the unstable length exceeds about 12 in. in aluminium 
alloy materials. For short cracks in very wide sheets, values of k calculated by 
equation (2) must still be valid. 
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3.3. SOLUTION OF BASIC EQUATION 


It is clear from equation (1) that, in order to derive a relationship between 
stress and crack length, knowledge of the material stress-strain curve is required. 
If the complete curve is available, having found the appropriate value of k, 
equation (1) can be solved for < and the corresponding stress f read from the curve. 

If the curve is not available, as is generally the case, and only normal material 
properties are known, use can be made of the generalised stress-strain relationship 
given in R.Ae.S. Structures Data Sheet 00.02.04. 
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+30 
26 
4-124 
22 | 
| +18 
be 
\ 
i4 
12 
+10 
+—+10 
FiGurE 5. Nomogram solution of equation (4). 
This relationship is 
m 
t 
Substitution into equation (1) with some transposition gives 
(4) i (1+ max ). @ 
fi k fe, 


Figure 5 is a nomogram solution of this equation. When values of f,, max, E and m 
are known for the material and k has been found from the appropriate figure, the 
net area stress f can be determined from the nomogram. This figure is reproduced 
from R.Ae.S. Structures Data Sheet 00.02.05 with the notation suitably modified. 
Reference should be made to this sheet and to Data Sheet 00.02.04 for detailed 
information on this method of solution. 


3.4. CORRELATION WITH EXPERIMENTAL RESULTS 


Figures 6 to 11 show curves of gross area stress against unstable crack length 
for several materials and sheet sizes. All curves were derived using values of k 
obtained from Figs. 2, 3 and 4 and stress-strain properties from control tests on 
the material from which the experimental results were obtained. Within the limita- 
tions of the available experimental data these figures show materials of markedly 
different characteristics. It is important to realise that differences in the relation- 
ship between stress and unstable crack length for different alloys of the same basic 
group result entirely from differences in their stress-strain curves. The effect of the 
stress-strain curve parameters and sheet size is illustrated in the Appendix by a 
detailed analysis of some experimental data. 
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Figure 10. Theoretical curve and experimental results for $.520 steel (b= 10 in.). 
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4. Discussion 


In spite of the approximate nature of the analysis it is seen to agree remarkably 
well with experimental data from a range of typical aircraft structural materials. 
This shows that the relationship between unstable crack length and stress in sheet 
materials depends on three factors. 


(1) The inherent “sharpness” of the fatigue crack—essentially a 
property of the material 


(2) The geometry of the sheet as it affects the strain distribution 


(3) The material stress-strain curve properties, principally elongation 
and proof/ ultimate ratio. 


The basic principle may be applied to more complex two-dimensional stress 
systems, such as that in a sheet structure with discrete reinforcement, the only 
fundamental requirement being sufficient knowledge of the stress system to estimate 
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Figure 12. Comparison of three aluminium alloys. Theoretical curves and experimental 
results for 10 in. wide sheets. 


k values. References 8 and 9 give this information for reinforced sheet structure. 
For thick plates and solid sections the influence of biaxiality on the failure condition 
must be considered”. 

While it is realised that this analysis is in no sense an advance on the more 
usual Griffith-Irwin analysis, it affords a much clearer picture of the parameters 
which influence the unstable crack propagation characteristics of materials. As 
indicated earlier, the introduction of an energy condition at instability is 
unnecessary if no preceding stable crack growth is considered. The present analysis 
has an additional slight advantage in being applicable over the whole stress range 
where the usual method is restricted to net area stresses below about I:I times 
the yield stress. 
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FIGURE 13. Comparison of three aluminium alloys. Theoretical curves and experimental 
results for 30 in. wide sheets. 


Appendix 


INFLUENCE OF SHEET SIZE AND STRESS-STRAIN CURVE PROPERTIES ON UNSTABLE CRACK 
LENGTH BEHAVIOUR 


The effect of sheet size and stress-strain curve properties can be shown by detailed 
analysis of the data given in Ref. 7 for the alloys 2024-T3(A), 2024-T3(B) and 2024-T36. 
The data were obtained from tests on sheets 30 100-064 in. and 90x 300-064 in 
Both 2024-T3 materials were within specification but were obtained from different 
suppliers. 

The relevant mean material properties obtained from three control tests on each 
material are quoted in Table I. In Table II the curve of theoretical stress against crack 
length is calculated for 10 in. wide sheets of each material and in Table III for 30 in. wide 
sheets. In these tables k values are obtained from Fig. 2 and the stress-strain relation 
from the mean of control test curves. Conversion from net area stress f to gross area 
stress f, is accomplished in plotting Figs. 12 and 13 which give the theoretical curves 
and experimental results for 10 in. and 30 in. wide sheets respectively. (A straight line 


November 1961 405 


Material 


E (Ib. /in.? x 10~-®) 


E. 


2024-T3 (A) 


WELBOURNE 


TABLE I* 
STRESS-STRAIN CURVE PROPERTIES 


2024-T3 (B) 


2024-T36 


10°3 

71:7 
0-14 
0°1470 
0°83 


WIDE SHEETS 


2024-T 36 
f 
0:00742 594 
000555 $554 0-771 
000437 44:7 0-624 
000400 0°573 
000389 400 0:°558 


2024-T 36 
f fit 
000439 447 
0-00325 336 0°468 
000257 268 0374 
000232 246 0343 
0:00226 23:7 0330 


10°6 10:2 
68:2 
0:19 0:16 
E+ 01964 01669 
olf 0-75 0:79 
TABLE II* 
RELATION BETWEEN STRESS AND CRACK LENGTH FOR 10 IN. 
Material 2024-T3 (A) 2024-T3 (B) 
x,,/b k f f f lh, 
0:05 19°8 0:00992 52:8 0:774 0:00843 S37 0:797 
01 26:5 000742 51:3 0741 000630 543 0:776 
0-2 336 000585 502 0:725 000497 492 0:704 
03 0:00536 492 0-711 000455 450 0644 
0:4 378 000520 488 0-705  0:00441 43:8 0-626 
TABLE III* 
RELATION BETWEEN STRESS AND CRACK LENGTH FOR 30 IN. WIDE SHEETS 
Material 2024-T3 (A) 2024-T3 (B) 
x,/6 k f tlh f 
0-05 33:5 000587 50:3 0-736 0:00498 49:2 0:704 
0-1 5-2 000434 43:7 0-641 0:00369 37:0 0°529 
0-2 572  0:00343 35:5 0°520; 0:00292 29:3 0-419 
03 63:2 0:00311 32:2 0-472 | 0:00264 266 0-381 
0-4 649 000303 0:00257 2$9 0:37) 
TABLE IV* 
RELATION BETWEEN STRESS AND CRACK LENGTH FOR 2 IN. WIDE SPECIMENS 
Material 2024-T3 (A) 2024-T3 (B) 
x,/b k f f lf, e f flf, 
6:05 9-5 | 0:02067 56:0 0821) 001757 59:1 0-845 
0-1 12°4 | 0:01584 542 | 001346 0-825 
0-2 15-4 | 001275 534 0°783 0:01084 56°5 0-808 
0:3 168 001169 52:8 0774, 0:00993 562 0-804 
0-4 172 001142 0:00970 56:0 0801 


0:773 


*Values of f are in Ib./in.? x 10-%, 
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2024-T36 
fi, 
001547 640 0893 
001185 626 
000955 61:0 0-851 
000875 60:4 
000855 60:2 0-840 
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joining some value of f, on the vertical axis to x,/b=1-0 on the horizontal axis is a line 
of constant net area stress, f=f,,.) 


Considering Fig. 12, it can be seen that there are two quite distinct regions of the 
curves of stress against crack length. For short crack lengths and correspondingly low k 
values, general yielding occurs before failure in all three materials and their relative 
strengths are determined by their proof/ultimate ratios. For longer cracks, the value of 
k has increased such that the net area stresses are still below the yield region at failure: 
consequently the relative strengths of the materials are determined by their respective 
elongations. 


In Fig. 13 the basic result is the same except that the general yielding region is 
reduced in size. Even so, the reversal of order of the three highest experimental points 
with respect to the elastic region confirms its effect. It is apparent that the increase in 
k values by a factor of approximately ./3 caused by the 3:1 increase in sheet width 
keeps net area stress at failure substantially elastic except for cracks less than about 
14 in. long. 


It follows from the previous paragraph that, if much narrower specimens, say 2 in 
wide, were tested, k values would be sufficiently low that general yielding would occur 
before failure, whatever the crack length. Thus, the relative strengths of the three 
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Figure 14. Comparison of three aluminium alloys. Theoretical curves for 2 in. wide 
specimens. 
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materials would be determined entirely by their proof/ultimate ratios, so that their 
apparent order of merit for crack resistance would be completely reversed to that which 
is observed in specimens of more practical dimensions. In Table IV, the theoretical 
curve of stress against crack length is calculated for 2 in. wide specimens of each 
material, k values being computed as in Section 3.2. The curves are plotted in Fig. 14 
and show the anticipated effect. 


The general conclusion from this is that comparisons of the unstable crack 
propagation characteristics of materials are incomplete without reference to specimen 
geometry and size or the amount of yielding of the uncracked section. The simple 
analysis proposed in this note shows the influence of these factors for the particular case 
of an unreinforced thin sheet. 
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